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1.  Introduction 


The  standard  nonlinear  optimization  problem  is  to  minimize  a  function  f(X)  that  depends 
on  a  design  variable  vector  X.  This  function  is  minimized  in  a  feasible  region  defined 
by  inequality  and  equality  constraints  g(X)  and  h(X),  respectively.  Recently,  a  more 
efficient  approach  has  been  developed  using  intermediate  design  variables  Y  and  response 
quantities  q(Y)  that  may  be  better  approximated  than  the  actual  functions,  but  still  can  be 
easily  related  back  to  the  original  functions  and  design  variables.  (Schmit,1977:19-24), 
(Barthelemy,1991),  and  (Starnes,  1979:564-570)  The  problem  now  is  to  minimize  the  new 
functions  related  io  f(X),  g(X},  h(X),  and  Q  based  on  the  better  approximated  design 
variable  and  response  functions.  To  carry  out  the  exact  calculation  of  the  response 
quantities  is  typically  quite  expensive,  so  they  are  normally  approximated  using  a  linear 
Taylor  series  approximation  in  the  intermediate  design  variable  Y  of  the  form: 

q(Y)  =  q(YQ)  +  Vq^iYQ)AY.  The  values  of  the  approximate  functions  are  then  calculated 
and  used  to  converge  by  iteration  to  a  suitable  tolerance  of  the  exact  solution.  Using  this 
iterative  approach  one  can  drastically  reduce  the  number  of  full  analyses,  provided  the 
approximations  are  of  high  enough  quality.  This  approach  and  many  similar  approaches 
all  use  the  linear  Taylor  series  approximation  to  the  response  function  based  on  a  single 
design  point.  New  approximations  are  constmcted  each  time  an  exact  analysis  is  made. 
Information  from  previous  analyses  is  typically  lost.  It  has  been  shown  that  higher 
accuracy  can  be  achieved  with  a  quadratic  Taylor  series  approximation  (Fleury,  1986:409- 
428);  however,  the  computational  expense  of  calculating  second  derivatives  often 
outweighs  the  benefit  gained  from  faster  convergence  of  the  more  accurate  estimates.  If  a 
computationally  less  expensive  method  could  be  developed  for  obtaining  approximate 
second  derivatives,  then  faster  convergence  should  be  obtained. 
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1.1  Objectives 

The  objective  of  this  thesis  is  to  test  a  method  of  numerical  optimization  that  uses  a 
quadratic  approximation  for  the  objective  function  and  constraints.  The  quadratic 
approximation  will  be  developed  using  a  Hessian  matrix  that  is  obtained  from  the  exact 
analyses  at  each  of  the  previous  points.  As  more  points  are  obtained  near  the  solution,  the 
approximation  for  the  Hessian  will  improve.  Also,  the  Hessian  will  be  calculated  weighing 
the  information  from  closer  points  more  heavily.  It  is  expected  that  using  this  quadratic 
approximation  will  improve  the  accuracy  of  the  approximations  that  are  sent  to  the 
optimizer,  resulting  in  faster  overall  convergence.  A  side  benefit  that  is  expected  from  the 
curvature  obtained  from  using  a  quadratic  approximation  is  that  the  dependence  of 
optimization  problems  on  move  limits  will  be  reduced.  A  second  quadratic  method 
Sequential  Quadratic  Approximation  (SQA)  will  also  be  tested. 


1.2  Problem  Statement 

The  overall  problem  is  to  find  a  way  to  mimic  a  second  order  Taylor  series  approximation 
without  actually  having  the  expense  of  calculating  second  derivatives.  The  method  used 
estimates  of  the  Hessian  matrix  by  using  information  that  is  already  available  from  previous 
analyses. 

1.3  Scope 

The  scope  of  this  thesis  is  to  develop  this  multipoint  quadratic  approximation,  then 
compare  it  to  other  methods  to  discover  if  it  can  solve  actual  problems  with  fewer 
iterations.  A  wide  range  of  stmctural  problems  have  been  chosen  for  comparison,  so  that 
some  idea  of  the  various  method’s  strengths  and  weaknesses  can  be  obtained.  These 
structural  problems  are  of  various  sizes  ranging  from  a  five  section  beam  with  one 
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constraint  to  a  40  member  frame  with  408  constraints.  These  problems  also  contain  a 
variety  of  constraints  including  stress  constraints,  displacement  constraints  and  frequency 
constraints. 
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2.  Background 

Typically  in  engineering  we  want  to  produce  the  best  design  possible  based  on  the 
constraints  given  to  us.  Optimization  concepts  have  provided  a  way  to  systematically 
produce  the  best  design  based  on  known  parameters.  A  typical  example  is  to  design  the 
lightest  structure  that  will  carry  a  prescribed  load.  The  weight  of  the  structure  is  the 
function  that  we  are  trying  to  minimize  and  the  allowable  stresses  and  displacements  are  the 
constraints.  We  can  make  the  structure  as  light  as  possible  but  we  must  do  so  in  a  manner 
that  does  not  cause  it  to  collapse.  In  the  past,  this  was  done  based  on  experience  and  trial 
and  error.  Numerical  optimization  techniques  have  now  given  engineers  systematic 
approaches  to  solving  these  problems. 

2.1  Optimization  Concepts 

The  best  way  to  explain  the  optimization  process  is  through  a  simple  example.  The  first 
step  in  the  optimization  problem  is  determining  what  is  to  be  minimized  (or  maximized); 
for  example,  minimizing  the  weight  of  a  structure.  In  order  for  the  problem  to  be 
optimized,  an  equation  or  function  must  be  developed  for  the  weight.  This  function  is  the 
objective  function,  which  will  be  called  f(X).  The  objective  function  is  a  function  of 
parameters  called  design  variables  X.  These  design  variables  typically  are  the  cross 
sectional  areas  of  the  members  of  a  structure  provided  that  the  lengths  of  the  members  are 
fixed.  The  next  step  is  to  develop  constraint  functions  also  based  on  the  same  design 
variables.  An  example  of  a  constraint  would  be  the  allowable  stress  in  a  member.  A 
function  must  be  developed  relating  the  cross  sectional  area  (the  design  variables)  to  the 
allowable  stresses.  These  are  our  constraint  functions.  There  are  also  constraints  called 
equality  constraints  that  relate  directly  to  the  design  variables.  Equality  constraints  are 
typically  explicit  functions  of  X  that  relate  one  of  the  design  variables  to  another.  An 
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example  of  the  would  be  a  constraint  where  the  cross  sectional  the  depth  of  a  beam  has  to 
be  twice  its  width.  Now,  we  have  an  objective  function /fXj,  constraint  functions  g{X),  and 
equality  constants  h(X)  with  X  being  the  vector  of  design  variables. 

The  objective  is  to  rniniinize/fX)  while  not  violating  any  of  the  constraints  in  g(X)  and 
h(X).  For  the  rest  of  this  thesis,  equahty  constraints  h{X)  are  not  used  in  the  numerical 
examples,  although  this  is  not  an  inherent  restriction.  The  next  question  is,  “how  do  we 
know  when  a  minimum  has  been  achieved”. 

The  three  necessary  and  sufficient  conditions  for  a  minimum  are  called  the  Kuhn-Tucker 
conditions.  (Vanderplats,1984;17)  These  are; 

1 .  X  is  feasible 

2. X.^g{Xy  =  Q  j=l,m  Aj>0 

m 

3.  V/(X)  +  £Vg^.(X)  =  0  Xj>0 

J=1 

If  a  constraint  is  not  active,  i.e.,  not  equal  zero,  then  its  A  is  zero.  Using  these 
conditions,  an  optimizer  solves  a  problem  given  a  function  and  it  constraints.  However,  it 
must  be  pointed  out  that  this  method  only  finds  a  local  minimum  of  a  function.  If  more 
than  one  local  minimum  exists,  the  optimizer  could  find  any  of  one  of  them.  There  is  no 
guarantee  that  the  one  it  finds  will  be  the  global  minimum.  For  this  thesis  it  is  assumed  that 
the  optimizer  can  do  its  job  effectively.  The  objective  in  this  thesis  is  to  determine  what  is 
the  best  form  to  express  the  objective  function  and  constraint  function  so  that  the  optimizer 
can  solve  the  overall  problem  quickly  and  efficiently. 


2-2 


2.2  Iterative  Procedure 


Figure  2-1:  Optimization  Using  Approximations 

For  structural  optimization,  the  overall  procedure  for  the  “approximation  concepts” 
approach  is  shown  in  Figure  2-1.  This  is  the  approach  used  by  many  commercial 
structural  analysis  programs  that  also  include  optimization,  such  as  MSC/NASTRAN, 
IDEAS,  and  ASTROS.  The  first  step  is  to  perform  the  initial  analysis  based  on  some 
starting  point.  It  is  generally  best  to  have  this  starting  point  be  a  feasible  solution  to  the 
problem  but  it  is  not  necessary.  Starting  at  a  highly  infeasible  point  far  from  the  optimum 
can  make  the  optimization  take  longer.  The  initial  analysis  must  give  us  the  value  of  the 
objective  function  and  the  values  of  all  the  constraints,  f(X)  and  g(X),  respectively.  The 
second  step  is  to  do  the  sensitivity  analysis  to  obtain  the  derivatives  of  the  objective 
function  and  constraints.  The  third  step  is  the  approximate  problem  generator,  the  focus  of 
this  thesis.  This  step  could  be  skipped  and  the  exact  problem  could  be  given  to  the 
optimizer  so  that  it  can  solve  the  problem  to  find  the  optimum.  The  disadvantage  with  this 
approach  is  that  this  typically  requires  many  analyses.  If  each  analysis  takes  a  lot  of 
computer  time  then  it  will  take  a  long  time  to  find  the  optimum.  A  better  way  is  to 
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approximate  the  exact  problem  with  an  explicit  equation  that  gives  us  nearly  the  same 
values  of  the  objective  function,  constraints  and  the  gradients  for  each.  These 
approximations  will  only  be  good  near  the  point,  therefore  move  limits  are  placed  on  the 
design  variables.  Move  limits  are  restrictions  on  how  much  the  optimizer  can  vary  the 
design  variables  to  find  the  optimum.  Move  limits  are  added  as  another  constraint.  For 
example,  if  we  give  the  optimizer  a  starting  point  of  two  design  variables  of  (1,2)  and 
impose  a  move  limit  of  20%,  this  means  that  the  optimizer  is  restricted  to  use  values  of  the 
design  variable  in  the  range  of  (.8-1.2,1.6-2.4).  Move  hmits  of  this  type  are  called  box 
move  limits.  Another  type  of  move  hmit  is  a  circular  or  spherical  move  limit  for  which 
changes  in  one  design  variable  limit  changes  in  another.  The  way  this  is  imposed  is  by 
taking  the  dot  product  of  the  normalized  design  variable  vector  with  itself  and  setting  this 
value  equal  to  some  parameter,  given  by 


AX^SAX  <  r 

where  AX  is  the  change  in  the  design  variable,  S  is  a  scaling  factor  matrix  to  normalize  the 
variables  and  r  is  the  move  limit,  such  as  0.2.  5  is  a  square  diagonal  matrix  with  one 
divided  by  the  square  root  of  the  original  x  value  on  its  diagonal  for  each  contraint.  For 
example,  if  X=(l,2)  then 


Axj 

Ax2 


\T 


0 

1 

V2, 


|L^2J 


<0.2 


From  this  one  can  see  the  size  of  the  change  in  one  variable  limits  the  change  in  the  size  of 
another.  Instead  of  being  a  box  this  is  the  equation  of  a  circle.  In  multidimensional  space 
the  box  becomes  a  cube  or  hypercube  and  the  circle  becomes  a  sphere  or  hypersphere.  One 
of  the  objectives  of  this  thesis  will  be  to  see  if  spherical  move  limits  have  an  advantage  over 
box  type  move  limits. 

Once  the  move  limits  are  established,  the  approximate  problem  is  then  given  to  the 
optimizer.  The  optimizer  then  proceeds  to  find  the  approximate  optimum  solution  within  the 
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bounds  of  the  move  limits  and  constraints.  Once  it  finds  this  intermediate  design  point,  the 
new  X  vector  is  used  in  the  next  step  where  a  detailed  analysis  is  performed  at  the  new 
point.  In  the  next  step  the  actual  objective  function  value  and  the  actual  constraints  are 
compared  to  the  approximate  optimum  received  from  the  optimizer.  If  these  are  within 
specified  tolerance  of  each  other,  then  the  optimization  process  is  complete.  If  not,  then  a 
new  sensitivity  is  computed  based  on  the  new  point  that  should  be  closer  to  the  optimum 
than  the  previous  point.  A  new  approximate  problem  is  then  created  and  this  is  given  to  the 
optimizer  to  solve.  Once  solved,  convergence  is  again  checked.  With  each  iteration,  the 
routine  should  get  closer  to  the  tme  optimum.  One  problem  often  encountered  is  that  the 
solution  gets  close  to  the  optimum  but  has  trouble  converging.  This  is  because  the 
accuracy  of  the  approximation  and  the  convergence  criteria  may  not  be  of  the  same  order. 
To  prevent  this  problem,  the  approximation  must  be  made  more  accurate.  This  can  be  done 
by  systematically  reducing  the  move  limit,  since  the  closer  we  are  to  the  point  of  the  last 
actual  calculation,  the  more  accurate  the  approximation.  The  key  factor  with  move  limits 
is  to  have  them  large  enough  when  optimization  begins  so  that  the  approximation  can  move 
closer  to  the  optimum  at  a  reasonable  rate,  and  having  them  small  enough  at  the  end  so  the 
approximation  is  accurate  enough  to  converge.  However,  as  the  accuracy  of  the 
approximation  increases  the  dependence  on  having  the  proper  (small  enough)  move  limits 
should  decrease. 

2.3  Assumptions 

This  thesis  focuses  strictly  on  step  three  in  Figure  2-1,  the  approximate  problem  generator. 
It  is  assumed  that  the  programs  that  calculate  the  analysis  and  sensitivities  are  correct. 
Also,  it  is  assumed  that  the  optimizer  does  its  job  efficiently.  The  optimizer  used  was 
“constr”  in  MATLAB  (Grace,  1993:9-11),  and  for  the  most  part  it  worked  well.  However, 
sometimes  it  had  trouble  solving  the  approximate  problem.  This  is  outside  the  scope  of  this 
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thesis  because  the  approximate  problems  developed  here  are  generic  and  could  easily  be 
adapted  to  work  with  any  such  optimizer. 
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3.  Literature  Survey 


Figure  2-1  shows  the  typical  optimization  process.  The  focus  of  this  thesis  is  on 
step  3,  which  is  generating  the  approximate  problem  for  the  optimizing  algorithm  to  solve. 
The  basic  concept  is  that  the  more  accurate  the  approximate  problem,  the  faster  the  problem 
will  converge  to  the  optimum. 

3.1  Early  Work 

In  the  1950’s  the  first  implementation  of  mathematical  programming  was  used  by 
Heyman  and  Livesly  to  optimize  a  single  bay  structure  subject  to  a  single  load  case 
(Kolonay, 1987:6).  It  quickly  became  obvious  that  optimizing  a  structure  involved 
numerous  iterations  and  would  involve  numerous  solutions  to  reanalyze  the  structure  and  to 
calculate  the  constraints.  In  the  1960’s,  Moses  developed  the  technique  of  using  a  Taylor 
series  expansion  to  approximate  the  objective  function  and  constraints.  (Kolonay, 1987:6) 
By  doing  this,  he  reduced  the  number  of  analyses  required  to  solve  the  optimization 
problem.  The  typical  optimization  problem  used  a  first  order  Taylor  series  approximation 
using  the  information  from  the  last  full  analysis.  A  first  order  approximation  was  used 
because  calculating  second  derivatives  needed  for  the  Hessian  matrix  for  the  second  order 
approximation  was  very  computationally  expensive. 

3.2  Recent  Work 

Recent  work  in  the  area  of  improving  the  approximation  used  by  the  optimizer  in 
stmctural  optimization  has  taken  several  directions.  The  first  direction  is  with  coming  up 
with  intermediate  design  variables  such  as  reciprocal  design  variables  that  can  form  better 
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approximate  nonlinear,  but  explicit,  functions  than  the  original  design  variables.  The 
second  direction  is  using  information  from  more  than  one  analysis  or  design  point.  An 
application  of  the  second  direction  is  to  use  previous  information  to  develop  a  quadratic 
approximation.  Pertinent  to  both  of  these  methods  is  the  question  of  move  limits.  Move 
limits  are  how  far  we  let  the  optimizer  go  in  using  the  approximation.  A  Taylor  series 
approximation  is  an  approximation  about  or  near  a  point;  therefore,  the  farther  away  from 
the  point  the  worse  the  approximation.  Determining  how  far  away  from  the  point  the 
approximation  is  adequate  is  a  science  and  art  all  to  itself.  However,  if  one  can  increase  the 
distance  where  the  approximation  is  valid,  this  should  also  increase  the  rate  of 
convergence. 

It  was  discovered  that  by  using  Taylor  series  approximations  in  reciprocal  variables,  the 
accuracy  of  the  approximations  improved  and  faster  convergence  was  achieved. 
(Schmit,  1974:692-699).  As  Vanderplaats  states,  “using  reciprocal  design  variables  can 
make  non-linear  structural  optimization  problems  linear,  so  that  a  Taylor  series 
approximation  is  very  accurate”  (Vanderplaats,1984,260-263).  This  technique  has  reduced 
the  number  of  analyses  required  in  many  structural  problems  and  it  worked  well  because 
the  quantities  being  reciprocated  (such  as  cross  sectional  areas)  are  positive  and  not  near 
zero.  In  some  optimization  problems  where  a  variable  could  be  zero,  the  reciprocal 
approximation  could  become  undefined. 

Typically  the  information  used  to  develop  the  approximation  was  only  based  on  the 
current  point  and  the  information  from  earlier  analyses  was  lost.  Haftka  and  coworkers 
initially  investigated  multipoint  approximations  with  mixed  results  (Haftka,1987:289-301). 
In  particular,  they  concluded  that  two  and  three  point  approximations  may  be  no  better  than 
the  reciprocal  approximations  for  structural  problems.  Fadel  incorporated  previous  gradient 
information  in  a  two  point  approximation  that  was  used  to  select  better  intermediate  design 
variables  for  a  linear  approximation  (Fadel,  1990, 117-124).  In  contrast,  Rasmussen 
incorporated  only  previous  function  values,  while  ignoring  gradients  from  previous 
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iterations  (Rasmussen,  1990:253-258).  He  demonstrated  stable  convergence  for  a 
cantilever  box  beam  where  other  methods  failed  to  converge.  Recently,  Wang  and  Grandhi 
demonstrated  more  rapid  convergence  for  frame  structures  by  employing  multivariate 
splines  based  on  aU  available  exact  analyses  (Wang,  1994: 2090-209 8).  However,  their 
method  does  not  always  reproduce  known  function  values.  In  another  paper  with  the  help 
of  Canfield  they  did  reproduce  function  and  gradient  values  (Wang,1994:269-280). 

Miura  and  Schmit  demonstrated  the  higher  accuracy  of  quadratic  Taylor  series 
approximations  for  stmctural  optimization  with  frequency  constraints  (Miura,  1978:336- 
351).  However,  the  computational  savings  due  to  faster  convergence  were  balanced  by  the 
increased  cost  of  calculating  second  derivatives.  Snyman  and  Stander  had  success  in  using 
a  two  point  successive  approximation  which  uses  a  pseudo  second  order  approximation 
based  on  these  two  points  (Snyman,  1994).  The  Hessian  calculated  for  this  case  was 
diagonal  with  aU  terms  on  the  diagonal  being  the  same  value.  The  intent  was  not  to 
estimate  the  exact  Hessian  but  to  give  the  optimizer  some  curvature  information  at  the  new 
point.  The  benefits  of  this  approach  were  reduced  dependence  on  move  hmits  and  slightly 
faster  convergence.  Also,  a  large  amount  of  Hessian  information  did  not  need  to  be  stored. 
In  1993  Toropov,  Filatov,  and  Polynkin  developed  a  method  of  using  multiple  point 
approximation  in  a  subregion  to  where  finite  element  analyses  are  performed,  then  exphcit 
approximations  of  the  objective  and  constraint  function  are  obtained  and  sent  to  the 
optimizer  to  be  optimized.  The  subregions  are  selected  based  on  the  move  limits.  Their 
method  uses  weighting  factors  that  give  values  closer  to  known  points  more  weight 
(T  oropov,  1993:7-14). 

3.3  Conclusion 

All  of  the  methods  discussed  above  are  a  step  towards  a  general  method  to  develop  a 
better  approximation  for  the  optimizer  to  use.  An  approach  that  could  combine  the  best  of 
these  could  possibly  be  more  general  and  robust.  If  one  could  use  reciprocal  design 
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variables  and  put  the  information  obtained  from  all  previous  analyses  to  use,  it  would  seem 
feasible  that  a  better  approximation  could  be  found  and  used  for  faster  convergence  of  the 
optimization  problem. 
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4.  New  Theory 


4.1  Overview 

The  focus  of  this  thesis  is  to  develop  and  test  a  new  approximation  method  to  be  used 
by  an  optimization  routine.  The  basic  concept  is  to  make  use  of  the  information  obtained  in 
each  analysis  to  provide  a  more  accurate  approximation  of  a  function.  Typical  optimization 
routines  use  a  first  order  Taylor  series  approximation  because  it  is  too  computationally 
expensive  to  use  the  actual  function.  Also,  the  approximation  typically  is  based  only  on  the 
latest  point.  The  information  from  previous  analyses  is  no  longer  of  any  use.  The 
approach  used  here  is  to  use  the  information  obtained  in  aU  previous  analyses  to  create  a 
quadratic  approximation  of  the  function.  This  is  done  without  calculating  second 
derivatives.  A  side  development  during  this  research  was  the  creation  of  spherical  move 
limits  that  could  be  used  with  any  approximation  routine. 

To  compare  whether  any  benefit  is  obtained  from  this  quadratic  approach,  it  will  be 
compared  to  other  methods  already  being  used.  All  of  the  optimization  routines  use  one  of 
four  types  of  approximations.  The  first  type  is  using  direct  variables  and  a  linear 
approximation,  the  second  is  using  direct  variables  and  a  quadratic  approximation,  the  third 
is  using  reciprocal  variables  with  a  hnear  approximation  and  the  fourth  is  reciprocal 
variables  with  a  quadratic  approximation. 

4.2  Direct  vs.  Reciprocal  Variables 

Direct  variables  are  usually  physical  quantities  in  the  optimization  problem,  such  as  the 
cross  sectional  areas  of  a  structural  member.  A  common  objective  function  in  a  structural 
problem  is  the  weight  of  a  structure  which  is  linear  in  the  cross  sectional  area.  This  is  the 
case  for  all  of  the  problems  in  this  research.  Since  the  objective  functions  are  linear,  a  first 
order  Taylor  series  approximation  is  exact.  However,  the  constraints  for  structural 
problems  are  usually  not  linear  functions.  A  typical  constraint  for  a  structural  problem  is  a 
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mayimnm  allowable  stress,  which  is  a  function  of  force  divided  by  area.  Therefore,  a  linear 
approximation  of  the  constraint  is  probably  not  very  accurate.  If,  instead  of  using  the  cross 
sectional  area  as  the  design  variable,  one  uses  an  intermediate  design  variable  that  is  equal 
to  one  divided  by  the  cross  sectional  area,  then  a  better  approximation  can  be  achieved.  The 
stress  then  becomes  a  linear  function  of  this  new  design  variable.  Such  variables  are  called 
reciprocal  variables  (Vanderplats,  1984:260-263).  Their  use  can  greatly  improve  the 
accuracy  of  a  Taylor  series  approximation  for  many  stmctural  problems.  However,  if  the 
values  of  the  design  variables  get  very  small  or  pass  through  zero,  using  reciprocal 
variables  can  result  in  functions  becoming  undefined.  In  these  cases,  reciprocal  variables 
will  not  be  very  useful. 

4.2.1  Direct  Variable  Approximations 

The  first  order  Taylor  series  approximation  of  the  multivariate  function /fx)  near  some 
point  Xq  is 

7^(x)  =  /(Xo)  +  [Vf(Xo)f  Ax  (4-1) 

with  X  being  a  vector  of  the  design  variables.  Two  of  the  methods  used  in  this  thesis  use 
this  approximation.  They  are  SPA-DS,  which  stands  for  single  point  approximation  using 
direct  variables  with  spherical  move  limits  and  SPA-DB  which  is  a  single  point 
approximation  using  direct  variables  with  box  move  limits.  SPA-DB  is  the  same  as  what  is 
usually  called  sequential  linear  programming  (SLP).  The  quadratic  Taylor  series 
approximation  contains  all  the  information  of  the  linear  plus  a  second  order  term.  It  is  of 
the  form: 

7e(x)  =  /(Xq)  +  [Vf(Xo)f  Ax  +  ^  Ax^HAx  (4-2) 

where  H  is  the  Hessian  matrix.  The  approach  tested  in  this  research  is  to  develop  a 
reasonable  approximation  for  the  tme  Hessian  based  on  previous  points.  One  method 
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tested  uses  this  approximation;  it  is  caUed  MQA-D,  which  stands  for  multipoint  quadratic 
approximations  using  direct  variables. 


4.2.2  Reciprocal  Variables 


To  develop  the  first  order  single  point  approximation  in  reciprocal  variables  we  begin 
with  the  first  order  Taylor  series  approximation  in  the  intermediate  variable  y 

/L(y)  =  /(yo)  +  [Vf(yo)fAy  (4-3) 

which  is  the  same  as  eq.  (4-1)  except  using  y  instead  of  x. 

For  reciprocal  variables  y  =  y(x),  where  y.  =  — ,  we  define  the  single  point 

X. 

approximation  with  respect  to  x,.  as 

//?(^)  =  fiiy^x))  (4-4) 

where 

fy(yo)  =  f(Xo)  (4-5) 

and 


i=l  ,=1  <^0/  <^0/ 


1  1 


yXj  Xq- j 


for  which  n  is  the  number  of  design  variables  (length  of  the  x  vector). 


^Oi  _ _ ^  _  ^2 

.2  ~  Xq,- 


^oi  yli 

Therefore  the  final  equations  for  the  reciprocal  linear  approximation  for  x  is 

'  1  1 

- .  (4-6) 

This  approximation  is  used  by  SPA-RS  and  SPA-RB  which  are  single  point 
approximations  using  reciprocal  variables  and  spherical  and  box  move  limits,  respectively. 


//eW  =  /(^o)  +  X^(-4) 

(=1  '^Oi 
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The  quadratic  approximation  for  reciprocal  variables  starts  with  the  linear  approximation 
and  adds  another  term: 


/gfi  (^)  =  /(^o )  +  S  ^  (-4  )j 

i=l  0^0/ 


1  1 


VX,  Xo,.y 


+  |a/H,A3.. 


(4-7) 


The  Hessian  approximations  can  estimate  the  Hessian  with  respect  to  y,  therefore  no  chain 
rule  is  required  when  converting  the  second  order  term  to  use  variables  in  x.  The  final 
equation  is 


fQR(x)  =  /(xq)  +  S^(-4)| 
,=1  COCoi 


0 

1 

0 

- — - 

-I-- 

- — - 

— 
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^X,.  XoJ 

yu 

^X,.  XoJ 

.  (4-8) 


Eq.  (4-8)  is  the  equation  used  by  MQA-R,  multipoint  quadratic  approximation  using 
reciprocal  variables  and  by  SQA,  the  sequential  quadratic  approximation. 


4.2.3  Using  The  Approximations 


Now  that  these  approximations  have  been  developed,  the  next  step  is  to  use  them  to 
solve  an  optimization  problem.  The  approximation  for  the  objective  function /(X)  is  f{X) 
and  could  be  any  of  the  four  previously  mentioned,  although  typically  the  approximation  of 
the  objective  function  is  linear  in  direct  variables.  The  approximation  for  the  constraint 
function  g(X)  is  g(X).  There  will  be  an  approximation  for  each  constraint.  The  problem 
now  becomes  one  of  finding 

/(x*)  =  min/(x),  xgQ  (4-9) 

^  =  (4-10) 

where  Q  is  a  feasible  region  where  none  of  the  constraints  are  violated.  These  equations 
are  solved  iteratively  until  the  approximations  converge  to  the  exact  solutions  to  within  a 
certain  tolerance. 
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4.3  Estimating  the  Hessian 


4.3.1  Using  All  Previous  Information 

Canfield  developed  the  approach  which  is  to  be  tested  (Canfield,  1994).  The  following 
parallels  his  derivation.  We  begin  with  a  quadratic  Taylor  series  approximation  of  a 
function  of  x. 

/g(x,H)  =  /(x^)  +  [Vf(Xj,)f  Ax  +  |ax^HAx  (4-11) 

Where  '>^k  is  the  design  vector  at  the  kth  iteration,  Ax=  x^^  -  is  the  change  in  the 
design  vector,  Vf(x^)  is  the  function  gradient,  and  H  is  an  approximation  to  its  Hessian. 
If  H  were  composed  of  analytic  second  derivatives,  eq.(4-ll)  would  represent  a  second 
order  Taylor  series  approximation.  Instead  of  calculating  the  actual  Hessian  we  will  use  the 
information  obtained  at  previous  points  to  calculate  the  Hessian  required  for  the  function  to 
pass  through  these  points  and  to  match  the  gradients  as  closely  as  possible.  After  each 
iteration  more  information  is  available  and  a  new  Hessian  is  calculated  using  the  additional 
information.  Also,  points  that  are  closer  to  the  current  point  are  given  more  weight  because 
closer  information  should  give  the  best  approximation.  As  we  get  closer  to  the  optimum 
the  approximate  Hessian  should  become  closer  to  the  actual  Hessian  and  should  speed 
convergence.  The  first  step  is  to  develop  a  set  of  equations  that  can  be  solved  to  find  the 
elements  of  the  Hessian.  We  will  begin  with  a  first  order  Taylor  series  approximation  as 
shown  eq.(4-12). 

=  /(Xfc)  +  [Vf(x,)f  (Xp  -  X,)  (4-12) 

p  e  P  =  l,2,...,k-1 

Usually  the  first  order  approximation  is  not  equal  to  the  actual  function  value.  However,  if 
these  two  values  are  known,  such  as  at  previous  points,  then  the  elements  of  the  Hessian 
can  be  considered  unknown  parameters  which  can  be  found  by  the  known  error.  Eq. 
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(4-13)  sets  the  quadratic  approximation  equal  to  the  known  function  at  the  previous  points. 
The  only  unknowns  are  the  elements  of  the  Hessian  matrix.  It  is  obvious  that  this  problem 
is  underdetermined  until  sufficient  points  are  obtained. 

/g(x^,H)  =  /(xp  =  /(xj  +  [Vf(x,)f  (Xp  -  X,)  +  ^(Xp  -  x,)''H(Xp  -  x,)  (4-13) 

peP  =  1,2,...,^:-! 

Considering  that  the  symmetric  Hessian  matrix  has  N=n(n+l)/2  unknown  elements  for  n 
design  variables,  we  need  more  equations  to  solve  for  these  unknowns.  Therefore,  we 
would  also  like  to  reproduce  the  gradients  from  previous  design  points  as  closely  as 
possible.  The  gradient  of  the  quadratic  approximation,  eq.  (4-14)  evaluated  at  the  previous 
design  points  is 


VJq  (Xp ,  H)  =  Vf(x^ )  -H  H(Xp  -  Xt)  =  Vf(Xp ),  /?  G  P . 


(4-14) 


Equations  (4-13)  and  (4-14)  constitute  a  set  of  linear  equations  in  the  N  unknown 
elements  of  the  Hessian  matrix  H.  However,  the  equations  are  in  general  inconsistent.  The 
equations  are  solved  using  the  Moore-Penrose  pseudo-inverse  which  provides  a  unique 
solution  with  the  smallest  norm  of  the  terms  in  the  Hessian. 

The  first  step  in  solving  (4-13)  and  (4-14)  is  to  write  them  in  standard  form  for  solving  a 
set  of  linear  equations.  Eq.  (4-14)  becomes 

Xph«  =  7,,  P^P  (4-15) 

where 


AxJ  0 
0  AxJ 


0  0 


0 

0 

AxJ 


(4-16) 
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(4-17) 


Xpis  an  nxn^  matrix  and 


^p=Xp-X, 

is  the  change  in  x  from  each  previous  point  to  the  current  point  where 


h'  =[[//„  (4-18) 

is  a  vector  of  the  rows  of  H.  The  right  hand  side  is  the  vector  difference  between  the  actual 
and  approximate  gradient  for  each  previous  design  vector: 


Yp=fiXp)-fL{Xp) 

The  next  step  is  to  put  eq.(4-13)  into  standard  form.  It  can  be  rewritten  as 

|Ax^X,h''=«,,  psP 


(4-19) 


(4-20) 


where  the  right  hand  side  is  equal  to  the  difference  between  the  acmal  function  value  and  its 
linear  approximation. 


5p^f(Xp)-f,ix^)  (4-21) 

Once  eqs.(4-13)  and  (4-14)  have  been  put  into  this  form,  they  are  combined  into  a  set  of 
linear  equations  of  the  form 

Xh''  =  b.  (4-22) 

where 


■X,' 

X, 


X 


p 


IX  \ 


(4-23) 


which  is  a  matrix  created  from  the  left  hand  sides  of  eqs.  (4-15)  and  (4-20)  and  for  which 


'AxfXi 

1  AX2X2 

2  : 
.Ax^Xp 


(4-24) 
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The  combined  right  hand  side  of  eq.  (4-22)  is 

b=Lr.’'  rl  -y,  ^  ('*-25) 

We  impose  the  condition  that  H  be  symmetric  through  the  linear  transformation: 

h''=Th^  (4-26) 

where  T  is  an  n^  x  N  linking  matrix  of  zeros  and  ones,  with  one  non-zero  entry  per  row 

and  where  is  an  N  x  1  vector  of  the  symmetric  terms  in  H. 

The  next  step  is  to  calculate  the  weighting  factors  so  that  closer  points  are  given  more 
weight  and  we  emphasize  matching  function  values  over  gradients.  D  is  the  matrix 
containing  all  the  weighting  factors. 

D  =  diag(W,  Wj  •••  W)  (4-27) 

The  Wp’s  are  the  weighting  factors  for  each  previous  point  and  W  is  the  weighting  for 
functions  over  gradients.  Thus 

Wp=^pl„x«  (4-28) 

in  which  I  is  a  nxn  identity  matrix  and  is  a  scalar  weight  inversely  proportional  to  the 
square  of  the  distance  from  the  current  point. 

/),=|Kf  (4-29) 

The  weighting  factor  for  function  values  is 

W  =  wdiag()3,  P2  "■  Pp)  (4-30) 


where 


(4-31) 
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Once  all  the  factors  are  caleulated  the  final  expression  for  the  approximate  Hessian  becomes 

h''=T(DXT)^Db  (4-32) 

The  symbol  t  is  used  to  represent  the  Moore-Penrose  pseudo  inverse. 

To  help  illustrate  the  process  of  estimating  the  Hessian,  we  will  estimate  the  Hessian  for 
a  simple  function  by  using  the  equations  and  discussion  above.  Begin  by  picking  a 
function,  three  previous  points,  and  one  current  point. 

The  function  is 


3x 


2 


The  current  point  Xj.is  (1,1).  The  three  previous  points  are: 


^.5^ 

^75^ 

1  1 

a.5" 

^p2  = 

,.75, 

II 

,1.5, 

The  function  values  at  these  points  are: 


/(;c,)  = -0.02265 
/(x^,)  = -1.0453 
f{Xp2)  = -036353 
/(Xp3)  =  0.31823 


The  gradient  of  f(X)  is: 
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v/(X)  = 


2\  —  +  x. 


The  gradient  at  each  point  is 


/"-0.257^  /"-1.03 

1.28  }  '^^<^''>=  5.12 


-1.03>  f-0.458^  f-0.114A 

5.12  J’  2  28  }  ^-^^^^^^“10.569  / 


From  the  above  function  and  gradients  an  approximation  to  the  Hessian  at  point  will 
be  made  using  the  following  equations,  where  P  is  the  number  of  previous  points,  which 
is  3  for  this  case.  Now,  we  will  put  eq.  (4-14)  in  standard  form 

xX  =  rp,  P^P  (4-33) 

using  eqs.  (4-16),  (4-11),  (4-18)  and  (4-19). 

From  eq.  (4-16) 


■-.5 

-.5  0  0  -.25 

-.25 

0 

0 

5  .5 

0  o' 

_  0 

0  -.5  -.5  0 

0 

-.25 

.  X„3  = 

-.25 

0  0 

.5  .5_ 

and  eq.(4-19)  we  get  these  quantities; 


rpi  = 


-0.772"\ 
3.84  J’ 


7p2  = 


-0.200^ 
0.996  J’ 


7.3  = 


0.143 

-0.711 


The  next  step  is  to  put  eq.  (4-13)  into  the  form: 

-Ax^XX  =  ^.  P^P 

2  p  p  P’  ^ 

using  eq.  (4-21).  Therefore: 

Ux,,)  =  -0.534,  Ux,2)  =  -0.278,  Ux,^)  =  0-489 
5^1  =-0.511,  5^2  =-0.0852,  5^3  =-0.170 
Now  the  entire  set  of  equations  can  be  put  in  the  form 


(4-34) 


4-10 


Xh''  =  b 


(4-35) 


where 


M25 

.125 

.125 

.125 

.0313 

.0313 

.0313 

.0313 

.125 

.125 

.125 

.125 

■  -.5 

-.5 

0 

0 

0 

0 

-.5 

-.5 

-.25 

-.25 

0 

0 

0 

0 

-.25 

-.25 

.5 

.5 

0 

0 

0 

0 

.5 

.5 

.125 

.125 

.125 

.125 

.0313 

.0313 

.0313 

.0313 

.125 

.125 

.125 

.125 

b 


f  -0.772 " 
3.84 
-0.200 
0.996 
0.143 
-0.711 


-0.511 


-0.0852 


1,-0.170  J 


By  imposing  the  condition  that  H  is  symmetric,  we  will  reduce  the  number  of  unknowns  to 
just  the  symmetric  parts  using  the  linking  matrix  T,  so  that  h*^  =  Th'^. 


1  0  0 
0  1  0 
0  1  0 
0  0  1 


The  next  step  is  to  calculate  the  weighting  factors  using  eqs.  (4-27)-(4-31). 
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From  eq.  (4-29) 


P„=2.  ^,,  =  8,  /i„=2 
These  are  then  scaled  so  that  the  smallest  is  equal  to  1,  Thus 

P„=i.  /5,2=4,  P,,=l 

Next,  the  weighting  factor  for  each  previous  point  is  calculated  using  eq.  (4-28). 


'1 

o' 

'4 

O' 

'1  0' 

0 

1_ 

II 

0 

4_ 

II 

0  1_ 

The  weighting  factors  to  weight  functions  more  than  gradients  are  calculated  from  eqs. 
(4-30)  and  (4-31).  The  intermediate  calculations  for  eq.  (4-31)  involve 


= 


= 


'1 

01 

■-0.772' 

■-0.772' 

0 

ij 

.  3.84  _ 

.  3.84  _ 

'4 

0' 

r-0.200' 

■-0.800' 

_0 

4 

[  0.996 

3.98  _ 

ri  01 

'0.143' 

'0.143' 

o 

0.711 

0.711 

leading  to  the  weighting  factor 


3.98 

.0852 


46.73 


The  terms  of  eq.  (4-30)  are 

=  46.73,  wPp2  =  186.96,  =  46.73 
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All  of  these  weighting  factors  are  put  together  in  a  diagonal  weighting  matrix  D. 


D  = 


1 

0 

0 

0 

0 

0 

0 

0 

0 


0 

1 

0 

0 

0 

0 

0 

0 

0 


0  0  0 

0  0  0 

4  0  0 

0  4  0 

0  0  1 

0  0  0 

0  0  0 

0  0  0 

0  0  0 


0  0  0  0 

0  0  0  0 

0  0  0  0 

0  0  0  0 

0  0  0  0 

10  0  0 
0  46.73  0  0 

0  0  186.96  0 

0  0  0  46.73 


We  now  have  all  the  pieces  to  calculate  the  approximate  Hessian  from  eq.(4-32). 


r  2.05 


-0.911 

-0.911 


-2.96 


2.05  -.911 

-.911  -2.96_ 


The  actual  Hessian  is 

■  1.03  -0.51' 

-0.51  -2.05  ■ 

Use  of  this  approximate  Hessian  in  a  quadratic  approximation  should  help  the  accuracy 
compared  to  a  linear  approximation,  because  the  terms  are  of  the  approximate  Hessian  are 
of  same  order  and  sign  as  the  true  Hessian. 

When  first  testing  the  multipoint  quadratic  approximation  using  this  method  to  obtain  the 
approximate  Hessian,  the  optimizer  would  make  large  steps  that  would  reduce  the  objective 
function  and  create  large  constraint  violations.  It  was  discovered  that  this  was  the  result  of 
the  approximate  Hessian  sometimes  having  large  negative  terms.  The  optimizer  would 
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proceed  in  the  direction  to  take  advantage  of  these  negative  terms.  This  would  cause  the 
second  order  term  to  become  a  large  negative  term  overwhelming  the  first  order  term.  To 
prevent  this  from  happening  a  limit  was  placed  on  the  size  of  a  negative  contribution  of  the 
quadratic  term  in  the  approximation.  A  limit  of  10  %  of  the  size  of  the  first  order  term  was 
used  and  this  greatly  improved  the  accuracy  of  the  approximation.  The  10%  limit  worked 
well  for  the  reciprocal  case  and  15%  was  used  for  the  direct  variable  case.  Further 
refinement  of  these  hmits  could  possibly  produce  further  improvements  in  accuracy  of  the 
approximation. 

4.3.2  Using  Two  Most  Recent  Points 

Calculating  the  Hessian  using  the  pseudo  inverse  method  and  all  previous  point  will  get 
quite  expensive  when  there  are  many  design  variables  and  constraints.  In  this  section 
another  method  developed  by  Canfield  will  be  summarized  very  briefly.  This  method 
should  prove  to  be  less  computationally  expensive  but  still  give  a  useful  approximation 
(Canfield,  1995). 

As  with  the  previous  method,  the  elements  of  the  function  vector  f(X)  will  be  estimated 
using  a  quadratic  Taylor  series  approximation  of  the  following  form: 

/g(x)  =  /(x^)  +  [Vf(x  Jf  Ax  +  i  Ax^HAx  (4-36) 

Where  is  the  design  vector  at  the  kth  iteration,  Ax=  x^  -  x^_j^  is  the  change  in  the 
design  vector,  Vf(Xk)  is  the  function  gradient,  and  H  is  an  approximation  to  its  Hessian. 
Again,  if  H  were  composed  of  analytic  second  derivatives,  this  would  represent  a  second 
order  Taylor  series  approximation. 

Consider  updating  the  Hessian  based  on  the  two  most  recent  design  points 

(4-37) 

where  is  chosen  to  satisfy  the  condition 
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-  W(Xk)  (4-38) 

One  solution  for  U^.  is  an  update  matrix  based  on  Powell’s  Modified  BFGS  rank  two 


update 


_  <1^ 
~  T 


where 


e  = 


q  s  s  H^s 


q  =  dy  +  (l-e)H,^s 


s  =  x,-x^_^ 


1  if  s^y  >  O.lH^s', 


(4-39) 


s^H,s-sy 


otherwise 


The  Hessian  is  further  updated  using  the  following  update  formula  to  match  the  previous 
function  value. 

H,  =  H^  +  U:  =  H,.,+U,  +  U:  (4-40) 

where 


u;=^ 


f(Xk-i )  -  )-{Wix^)Ys-^s^H*s 


(s^sf 


ss 


(4-41) 


Hk  is  the  new  Hessian  approximation. 

4.4  Spherical  Move  Limit 

It  is  convenient  to  apply  a  spherical  move  limit  because  it  is  also  a  quadratic  constraint. 
Numerical  experiments  here  have  shown  that  convergence  seems  to  be  less  sensitive  to 
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spherical  move  limits  than  the  conventional  infinity  norm  side  constraints  in  which  the 
approximation  is  limited  to  a  hypercubical  region.  The  spherical  move  limit  is  posed  as 

(x  -  Xf^fSi^ix - Xj) -  <  0  (4-42) 

where  the  terms  of  the  diagonal  scaling  matrix  are  given  by 


c(t)  _ 
^ij 


i  =  l,2,...,n 


(4-43) 


at  the  Ml  iteration  for  n  design  variables.  The  radius  of  the  sphere  may  vary  from  one 

iteration  to  the  next;  however,  the  use  of  quadratic  constraint  approximations  can  decrease 
the  need  to  gradually  reduce  it,  as  is  often  required  for  first  order  approximations  to 
converge. 

For  most  of  the  cases  the  initial  move  limit  is  applied  and  a  reduction  factor  is  used  after 
each  iteration  to  decrease  the  move  limit  as  the  approximation  gets  closer  to  the  optimum. 
This  helps  convergence  by  not  letting  the  approximation  go  so  far  that  inaccuracies  keep  it 
away  from  the  optimum.  This  reduction  factor  can  apply  to  any  type  of  move  limit, 
spherical  or  box. 
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5.  Numerical  Results 


5.1  Overview 

The  purpose  of  this  research  is  to  develop  a  more  efficient  method  of  structural 
optimization  which  will  gamer  the  higher  accuracy  and  quicker  convergence  of  a  quadratic 
Taylor  series  approximation  without  incurring  the  expense  of  calculating  the  Hessian 
matrix  The  approach  involves  fully  utilizing  the  information  from  previous  iterations.  AH 
the  exact  stmctural  and  sensitivity  analyses  will  be  used  to  estimate  the  Hessian  matrix  for  a 
quadratic  approximation  to  the  response  quantities. 

Once  this  method  was  programmed,  it  was  tested  using  several  known  structural 
examples  from  hterature:  a  five  section  beam,  a  four  member  frame,  a  ten  bar  truss,  a 
twelve  member  frame,  a  113  member  tmss  called  ACOSS  n  linked  to  6  design  variables 
and  31  design  variables,  and  finally  a  40  member  frame.  These  seven  structures  have  a 
wide  range  of  constraints  including  stress,  displacement  and  frequency  constraints.  These 
structures  were  optimized  using  eight  different  techniques. 

The  techniques  are  broken  into  two  categories:  direct  and  reciprocal  variables.  The  first 
four  techniques  used  direct  variables  which  are  the  cross  sectional  areas  of  the  members. 
The  first  of  the  direct  methods  used  is  the  optimization  routine  in  MATLAB  called 
“constr”,  written  for  constrained  optimization  (Grace,1993:9-ll).  It  is  given  the  exact 
problem  to  solve.  This  method  is  also  used  by  all  the  others  to  solve  the  approximate 
problems.  The  second  direct  method  used  was  MQA-D  which  stands  for  Multipoint 
Quadratic  Approximation  in  Direct  variables.  This  method  is  an  implementation  of  the 
equations  in  section  4.3.1  used  to  generate  the  approximate  Hessian.  The  third  method, 
SPA-DB,  is  a  single  point  linear  approximation,  identical  to  Sequential  Linear 
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Programming  (SLP).  The  fourth  direct  method  was  SPA-DS  which  is  the  same  as  SPA¬ 
DE  except  it  uses  spherical  move  limits. 

There  were  four  reciprocal  variable  methods  used  for  comparison.  The  first  is  MQA-R, 
which  is  similar  to  MQA-D  above  except  that  reciprocal  intermediate  design  variables  are 
used.  It  is  compared  to  three  other  methods;  the  first,  SPA-RS,  is  a  single  point 
approximation  in  reciprocal  variables,  with  spherical  move  limits  discussed  in  the  section 
4.3.  The  next  method,  SPA-RB,  is  the  same  as  SPA-RS  except  that  box  type  move  limits 
are  used.  The  last  method,  SQA-R,  is  sequential  quadratic  approximation  designed  to  be 
less  computationally  expensive  than  MQA-R.  It  is  an  implementation  of  the  equations  in 
section  4.3.2. 

The  strucmral  analyses  for  the  4-Member  Frame,  12-Member  Frame,  and  40  Member 
Frame  are  done  by  a  program  called  FRAMIN  (Kolonay,1987).  The  ten  bar  truss  and  both 
the  ACOSS  models  are  analyzed  by  a  program  called  SALVO  (Canfield,  1986).  The 
optimization  routine  is  written  in  MATLAB  and  interfaces  with  SALVO  and  FRAMIN. 

CPU  times  were  measured  for  some  of  the  problems.  It  must  be  empasized  that  this  was 
done  just  to  get  a  preliminary  idea  of  the  expense  of  the  different  approximations.  It  is 
understood  that  estimating  the  Hessian  using  the  pseudo  inverse  can  get  computationally 
expensive.  The  intent  was  to  get  an  idea  of  how  large  a  problem  would  make  the 
computational  expense  unbearable  compared  to  the  expense  of  the  finite  element  analysis. 
None  of  the  structures  optimized  in  the  research  are  computaionally  expensive  to  analyze. 
Therefore,  it  is  doubtful  that  the  quadratic  methods  would  require  less  CPU  time.  The 
purpose  is  to  see  if  the  quadratic  methods  can  optimize  the  stmctures  in  fewer  iterations  and 
to  see  if  the  CPU  times  are  reasonable.  If  the  quadratic  method  can  optimize  a  stmcture  in 
fewer  iterations  then  perhaps  it  will  be  more  economical  for  a  certain  range  of  problems 
with  a  small  number  of  design  variables,  but  that  require  a  large  amount  of  CPU  time  to 
analyze.  Another  reason  why  the  CPU  times  should  only  be  used  for  a  preliminary  look  is 
that  the  finite  element  analyses  are  all  in  FORTRAN  which  is  very  fast  compared  to 
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MATLAB  in  which  the  optimization  routines  were  written.  Again,  the  CPU  times  should 
only  be  compared  to  see  if  the  method  is  in  the  ballpark. 

5.2  Comparison  of  Methods 

5.2.1  Five  Section  Beam 


Figure  5-1:  5-Section  Beam 
5. 2. 1. 1  Description 

The  first  model  optimized  is  a  5  section  cantilever  beam  obtained  from  a  paper  by 
Toropov,  Filatov,  and  Polynkin;  see  Figure  5-1  (Toropov,1993,7-14).  The  beam  is  a 
square  box  section  with  five  different  cross  sections.  The  thickness  of  each  section  is 
uniform  while  the  design  variable  Xj  controls  the  height  and  width  of  the  jth  section.  The 
constraint  is  the  tip  displacement.  The  objective  function  f(X)  and  constraint  function  g(X) 
are: 
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A  feasible  starting  point  Xo={5,5,5,5,5}  was  chosen  where /(Xoj=  1.560  and  g(XJ=0.0. 
Toropov  et  al.  converged  to  the  optimized  objective  function  of  1.34  in  nine  iterations 
using  convergence  criteria  of  the  change  in  objective  function  of  <  0.01  and  maximum 
constraint  violation  of  0.02. 

5.2. 1.2  Results 

5. 2. 1.2.1  Optimization  Problem  Overview 

Table  5-1  gives  an  overview  of  the  information  used  in  this  optimization  problem.  All  of 
the  methods  converged  to  the  same  optimum. 

Table  5-1  5-Section  Beam  Optimization  Overview 


No.  of  Design  Variables 

5 

No.  of  constraints  (type) 

l(Disp) 

Convergence  Criteria 

Max.  Change  in  Objective 

.001 

Function 

Max.  Constraint  Violation 

.001 

Max.  Change  in  Design  Variable 

.01 

Move  Limits 

.5  with  .9  reduction  factor 

Min.  Allowable  Design  Variable 

.01 

Max.  Allowable  Design  Variable 

10 

Initial  Objective  Function 

1.56 

Initial  Design  Variables 

Initial  Max.  Constraint  Violation 

0 

Optimum  Objective  Function 

1.34 

Optimum  Design  Variable  Values 

6.05,5.27,4.47,3.53,2.15 
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5-SECTION  BEAM 


5-SECTION  BEAM 


Figure  5-2:  Direct  Method  Optimization 
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5.2. 1.2.2  Direct  Variable  Optimization 

Figure  5-2  shows  the  results  for  the  four  methods  using  direct  design  variables.  Table  5- 
2  shows  the  exact  number  of  iterations  each  method  took  to  reach  the  convergence  criteria. 
This  chart  also  shows  the  CPU  time  each  method  required  to  converge. 

Table  5-2:  Direct  Method  Convergence 


Method 

No.  Of  Iterations  to  Converge 

Time  In  CPU  Sec 

constr 

16 

4 

MQA-D 

11 

50 

SPA-DS 

20 

39 

SPA-DB 

14 

13 

Since  this  is  a  very  simple  problem,  solving  the  exact  problem  is  very  quick.  It  is  very 
reasonable  for  ‘constr’  to  solve  this  problem  very  quickly  in  CPU  seconds.  However,  if 
this  were  not  the  case  and  it  took  much  longer  to  solve  the  exact  function  and  gradients  then 
MQA-D  has  a  considerable  advantage.  MQA-D  is  intended  to  be  more  CPU  efficient  for 
problems  where  the  analysis  is  very  expensive  which  is  definitely  not  the  case  for  this 
problem.  The  purpose  here  is  to  show  that  it  can  solve  the  problem  in  fewer  iterations  so 
that  computational  savings  could  be  achieved  if  the  analysis  were  very  expensive.  All  of 
these  methods  converged  in  spite  of  oscillation  as  they  approached  the  optimum. 
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5-SECTION  BEAM 


CMCO'^t  lOCONCOOO 


NO.  OF  ITERATIONS 


5-SECTION  BEAM 


Figure  5-3:  Reciprocal  Method  Optimization 


5- 


CO 
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5.2. 1.2.3  Reciprocal  Variable  Optimization 


Table  5-3  shows  the  comparison  of  the  four  reciprocal  methods.  Table  5-3  gives  the 
exact  number  of  iterations  to  converge  to  the  required  tolerance  and  the  CPU  time. 

Table  5-3:  Reciprocal  Method  Convergence 


Method 

No.  Of  Iterations  to  Converge 

Time  In  CPU  Sec 

MQA-R 

6 

6 

SPA-RS 

16 

11 

SPA-RB 

12 

11 

SQA-R 

6 

4 

The  two  quadratic  methods  are  very  comparable  in  solving  the  problem  and  solve  it  very 
quickly.  SPA-RS  gets  very  close  to  the  answer  in  four  iterations  but  still  takes  another  12 
to  meet  the  criteria.  If  the  criteria  had  been  less  stringent  it  would  have  done  almost  as  well 
as  the  quadratic  methods.  The  apparent  benefit  of  the  quadratic  approximations  is  then- 
ability  to  converge  without  sensitivity  to  move  limits.  In  contrast,  the  single  point  methods 
cannot  converge  to  the  required  tolerance  until  the  move  limits  are  reduced  enough  to  be 
compatible  with  the  termination  criteria.  SPA-RB  Converged  faster  than  SPA-RS  but  it 
was  not  closer  to  the  optimum  until  10  iterations. 

5.2. 1.2.4  Comments 

For  this  problem  the  quadratic  methods  definitely  helped  speed  convergence  as  far  as 
number  of  iterations,  and  for  the  reciprocal  case  it  is  faster  even  for  this  small  problem. 
Even  though  the  box  move  limit  single  point  approximation,  SPA-RB  converged  in  fewer 
iterations  than  the  spherical  move  limit  SPA-RS,  it  took  a  long  route  to  get  there.  With 
minor  adjustments  in  the  move  limits  and/or  convergence  criteria  the  spherical  move  limits 
would  be  superior.  It  is  interesting  to  note  that  MQA-D  took  almost  five  times  as  long  to 
solve  the  problem  in  CPU  seconds  but  only  twice  the  number  of  iterations  as  MQA-R. 
This  is  because  as  more  points  are  accumulated  it  takes  each  iteration  longer  to  calculate  the 
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approximate  Hessian.  The  CPU  times  are  not  considered  significant  for  this  exphcit 
problem,  which  did  not  involve  any  finite  element  analysis.  Rather  CPU  times  were 
measured  to  verify  the  approach  for  a  problem  with  an  analytical  solution. 


5.2.2  Four  Member  Frame 


5. 2. 2.1  Description 

The  four-member,  single  bay,  single  story  frame  with  I  cross  section  in  Figure  5-4  was 
subjected  to  three  load  cases  as  shown  in  Table  5-4.  The  structure  has  5  nodes  and  was 
subject  to  45  displacement  constraints  and  12  stress  constraints  for  the  three  load  cases 
shown  in  Table  5-2.  For  this  structure  there  is  a  fixed  relationship  defined  between  the 
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moment  of  inertia(I),  section  modulus  (S)  and  the  cross  sectional  area  (A):  I=0.20720xA^ 
S=0.39300xAl  Where  the  constant,  0.20720,  has  units  of  in' ^  and  the  constant,  0.39300, 
has  units  of  m  \  This  model  was  obtained  from  Ray  M.  Kolonay,  Wright  Labs,  WPAFB, 
Ohio.  Table  5-5  and  Table  5-6  give  material  and  element  properties. 


Table  5-4:  4-Member  Load  Cases 


Direction  of 

oad 

Load  Condition 

Node 

X 

JL _ 

1 

2 

-10000 

3 

-10000 

4 

-10000 

2 

2 

5000 

-10000 

3 

-10000 

4 

-10000 

3 

2 

-10000 

3 

-10000 

4 

-5000 

-10000 

Table  5-5:  4-Member  Frame  Material  Properties 


Material 

Steel 

Modulus  of  Elasticity 

29x10®  psi 

Specific  Weight 

.283  Ib/in" 

Upper  Limit  on  Area 

88.0  in^ 

Lower  Limit  on  Area 

.5  in'^ 

Displacement  hmits 

x(.3  in0,y(.15in0,z(50.0rad) 

Number  of  Loading  Conditions 

3 

Table  5-6:  4-Member  Element  Properties 

Area 

Section  Modulus 

Moment  of  Inertia 

Maximum  Allowable  Stress 

Initial  Weight 

5016.0  lb 
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5.2.2. 1.1  Optimization  Problem  Overview 

Table  5-7  gives  an  overview  of  the  information  used  in  this  optimization  problem.  All  of 
the  methods  converged  to  the  same  optimum. 

Table  5-7:  4-Member  Optimization  Overview 


No.  of  Design  Variables 

4 

No.  of  constraints  (type) 

45(Disp),  12(Stress) 

Convergence  Criteria 

Max.  Change  in  Objective 

.1 

Function 

Max.  Constraint  Violation 

.001 

Max.  Change  in  Design  Variable 

.1 

Move  Limits 

1  with  .97  reduction 

Min.  Allowable  Design  Variable 

.5 

Max.  Allowable  Design  Variable 

88 

Initial  Objective  Function 

5016 

Initial  Max.  Constraint  Violation 

-.702 

Optimum  Objective  Function 

3204 

Optimum  Design  Variable  Values 

12.72,12.92,12.92,12.72 
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Figure  5-5;  Direct  Variable  Optimization 


4-MEMBER  FRAME 
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5. 2. 2. 1.2  Direct  Variable  Optimization 

Figure  5-5  shows  the  results  for  the  three  methods  using  direet  design  variables.  Table 
5-8  shows  the  exact  number  of  iterations  each  method  took  to  reach  the  convergence 
criteria. 

Table  5-8:  Direct  Method  Convergence 

Method 
constr 

MQA-D 

SPA-DS 

SPA-DB 

SPA-DS  was  very  close  to  the  optimum  at  116  iterations.  It  met  the  change  in  objective 
function  and  the  maximum  constraint  violation  criteria,  but  it  did  not  quite  meet  the  change 
in  the  design  variable  criteria  when  it  was  stopped.  The  quadratic  method  solved  the 
problem  in  the  fewest  number  of  iterations,  but  did  not  have  a  tremendous  advantage  over 
‘constr’. 
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Figure  5-6:  4-Member  Frame  Reciprocal  Method  Optimization 
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5. 2. 2. 1.3  Reciprocal  Variable  Optimization 

Figure  5-6  shows  the  comparison  of  the  four  reciprocal  methods.  Table  5-9  gives  the 
exact  number  of  iterations  to  converge  to  the  required  tolerance. 

Table  5-9:  4-Member  Frame  Convergence 


Method 

No.  Of  Iterations  to  Converge 

MQA-R 

5 

SPA-RS 

14 

SPA-RB 

12 

SQA-R 

12 

All  of  the  reciprocal  methods  came  very  close  to  the  optimum  in  about  4  or  5  iterations; 
however,  the  single  point  methods  took  longer  to  converge  within  tolerance.  The  evident 
advantage  of  the  MQA-R  was  its  capability  to  converge  cleanly  without  tightening  of  the 
move  limits,  whereas  the  other  methods  needed  help  from  the  tightened  move  limits  to 
converge.  For  this  problem  the  reciprocal  methods  showed  a  tremendous  advantage  over 
the  direct  methods. 
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5.2.3  Ten  Bar  Truss 


5. 2. 3.1  Description 

The  ten  bar  truss  shown  in  Figure  5-7  is  a  classic  example  used  to  test  optimization 
routines  (Haftka,  1992:238-239).  It  has  ten  members,  six  nodes,  and  the  allowable  stress 
is  25,000  psi  in  tension  or  compression.  The  area  of  each  member  was  initially  10.0  in^. 
At  node  2  and  4  there  is  a  100,000  lb  load  in  the  negative  y  direction.  The  structure  was 
restrained  in  the  z-direction. 


5-16 


5. 2. 3. 2  Results 


5. 2. 3. 2.1  Optimization  Problem  Overview 

This  problem  has  an  interesting  occurrence  in  that  there  are  two  distinct  local  optima 
very  close  to  each  other.  Any  nonlinear  problem  can  have  more  than  one  minimum.  For 
this  thesis  either  is  considered  correct,  because  all  of  the  optimization  techniques  available, 
only  have  the  capability  of  finding  a  local  minimum. 


Table  5-10:  Ten  Bar  Truss  Optimization  Overview 


No.  of  Design  Variables 

10 

No.  of  constraints  (type) 

14(Disp) 

Convergence  Criteria 

Max.  Change  in  Objective 

1 

Function 

Max.  Constraint  Violation 

.001 

Max.  Change  in  Design  Variable 

.1 

Move  Limits 

1.5  with  .9  reduction  for  direct 

2.0  with  no  reduction  for  recip 

Min.  Allowable  Design  Variable 

.1 

Max.  Allowable  Design  Variable 

40 

Initial  Objective  Function 

6294.7 

Initial  Design  Variables 

20,20,20,20 

Initial  Max.  Constraint  Violation 

.3132 

Optimum  Objective  Function 

5060.9,  5076.7  * 

Optimum  Design  Variable  Values 

[30.52,  .1,  22.20,  15.22,  .1, 

.56,  7.46,  21.03,  21.53,  .1] 

[30.73,.!, 23.94,14.73,.!, 
.1,8.54,20.95,20.84,.!]* 

*  Two  different  optima. 
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Figure  5-8:  Ten  Bar  Truss  Direct  Optimization 
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5. 2. 3. 2. 2  Direct  Variable  Optimization 

Figure  5-8  shows  the  results  for  the  three  methods  using  direct  design  variables.  Also, 
shown  in  Table  5-11  are  the  exact  number  of  iterations  each  method  took  to  reach  the 
convergence  criteria. 

Table  5-11:  Ten  Bar  Truss  Convergence 


Method 

No.  Of  Iterations  to  Converge 

constr 

34 

MQA-D 

35 

SPA-DS 

50-h 

SPA-DB 

73+ 

For  this  problem  ‘constr’  and  MQA-D  are  very  comparable  in  final  efficiency;  however, 
MQA-D  appears  to  be  very  close  to  a  reasonable  engineering  solution  in  about  20  iterations. 
SPA-DS  and  SPA-DB  have  a  great  deal  of  trouble  because  of  the  move  limits.  The  move 
Umit  reduction  factor  closes  in  on  SPA-DS  and  SPA-DB  before  they  are  very  close  to  the 
optimum,  so  they  get  bound  from  reaching  it.  This  problem  is  very  sensitive  to  the  move 
limits  in  the  direct  variables.  The  same  move  limit  as  the  reciprocal  approach  could  not  be 
used  because  none  of  the  direct  methods  would  even  come  close  to  converging.  SPA-DS 
will  converge  in  about  50  iterations  if  the  move  limits  are  reduced  at  a  slower  rate,  such  as 
0.95. 
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Figure  5-9:  Ten  Bar  Truss  Optimization 


5-20 


5. 2. 3. 2. 3  Reciprocal  Variable  Optimization 

Figure  5-9  shows  the  comparison  of  the  four  reciprocal  methods.  Table  5-12  gives  the 
exact  number  of  iterations  to  converge  to  the  required  tolerance. 

Table  5-12:  Ten  Bar  Truss  Convergence 


Method 

No.  Of  Iterations  to  Converge 

MQA-R 

14 

SPA-RS 

14 

SPA-RB 

14 

SQA-R 

36+ 

A  large  move  limit  was  used  for  this  problem  so  that  the  accuracy  of  the  reciprocal 
approaches  could  be  tested.  MQA-R,  SPA-RS,  and  SPA-RB  solve  the  problem  in  about 
ten  iterations  and  then  took  about  three  more  to  reach  the  criteria.  SQA-R  had  trouble  with 
this  problem.  It  would  eventually  converge  in  about  150  iterations.  The  problem  here 
seems  to  be  that  SQA-R  uses  a  Hessian  approximation  that  is  positive  definite.  This  causes 
this  approach  to  be  very  conservative  if  the  true  Hessian  is  not  close  to  positive  definite.  I 
believe  that  this  is  what  is  happening  here. 

The  reciprocal  approaches  for  this  problem  have  an  advantage  over  the  direct  in  that  a 
move  limit  is  almost  not  necessary  for  the  reciprocal.  Also  the  reciprocal  methods 
converged  much  faster  when  getting  close  to  the  optimum.  For  this  problem  the  quadratic 
terms  in  reciprocal  approximations  do  not  speed  convergence. 
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5.2.4  12  Member  Frame 


5. 2. 4.1  Description 

The  twelve  member,  two  story,  single  bay  frame  with  I-section  elements  in  Figure  5-10 
was  subjected  to  three  load  cases  shown  in  Table  5-13,  stress  constraints  on  each  element 
and  limited  nodal  displacement  in  the  x-direction.  There  are  12  nodes  in  the  structure. 
Nodes  11  and  12  are  restrained  from  displacement.  For  this  stmcture  there  is  a  fixed 
relationship  between  the  moment  of  inertia(I),  section  modulus  (S)  and  the  cross  sectional 
area  (A):  I=0.20720xA^  S=0.39300xA^  (Kolonay,1987,92),  where  the  constant, 

0.20720,  has  units  of  in'^  and  the  constant,  0.39300,  has  units  of  in  *.  Table  5-14  and 
Table  5-15  give  element  and  member  properties. 
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Table  5-13:  12  Member  Frame  Load  Cases 


Direction  of  load 


x(lb) 


-5000 


-10000 


-10000 


-10000 


-5000 


-10000 


-20000 


-20000 


-20000 


-10000 


-5000 


-10000 


-10000 


-10000 


-5000 


-10000 


-20000 


-20000 


-20000 


-10000 


-5000 


-10000 


-10000 


-10000 


-5000 


-10000 


-20000 


-20000 


-20000 


-10000 


Table  5-14:  12  Member  Frame  Material  Properties 


Material 


Modulus  of  Elasticit 


Specific  Weight 


Upper  Limit  on  Area 


Lower  Limit  on  Area 


Number  of  Loading  Conditions 


Displacement  limits 


Node  1,6 


Node  3,8 


Steel 


29x10°  psi 


.283  Ib/in^ 


88.0  in" 


.5  in" 


3 
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Table  5-15:  Element  Properties 


Area 


Section  Modulus 


Moment  of  Inertia 


Maximum  Allowable  Stress 


Initial  Weight 


129000  psi 


13244.4  lb 


Table  5-16  gives  an  overview  of  the  information  used  in  this  optimization  problem.  All 
of  the  methods  converged  to  the  same  optimum. 

Table  5-16:  12  Member  Frame  Optimization  Overview 


No.  of  Design  Variables 


No.  of  constraints  (type) 


Convergence  Criteria 


12(Disp),  36(Stress) 


Max.  Change  in  Objective  1 
Function 


Max.  Constraint  Violation 


Max.  Change  in  Design  Variable  I  .  1 


Move  Limits 


Min.  Allowable  Design  Vanable 


Max.  Allowable  Design  Variable 


Initial  Objective  Function 


Initial  Design  Variables 


Initial  Max.  Constraint  Violation 


Optimum  Objective  Function 


1  with  .9  reduction 


13244 


30,30,30,30,30,30,30,30,30,30,30,30 


Optimum  Design  Variable  Values  14.15,  14.00,  13.94,  14.18,  14.50, 

21.15,  16.72,16.50,21.21,16.18,16.03 
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Figure  5-11:  12  Member  Frame  Direct  Optimization 
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5.2.4. 1.1  Direct  Variable  Optimization 

Figure  5-11  shows  the  results  for  the  three  methods  using  direct  design  variables.  Also 
shown  in  Table  5-17  are  the  exact  number  of  iterations  each  method  took  to  reach  the 
convergence  criteria.  This  chart  also  shows  the  CPU  time  each  method  required  to 
converge. 

Table  5-17:  12  Member  Frame  Convergence 


Method 

No.  Of  Iterations  to  Converge 

Time  In  CPU  Sec 

constr 

34 

31 

MQA-D 

17 

392.4 

SPA-DS 

46 

228.4 

SPA-DB 

51 

94.7 

MQA-D  did  very  well  as  compared  to  the  other  direct  methods  based  on  iterations  in 
solving  this  problem;  however,  as  shown  in  the  table,  it  took  a  great  deal  more  CPU  time 
than  the  others.  The  analysis  is  done  so  quickly  that  the  cost  of  calculating  the  approximate 
Hessian  takes  up  the  majority  of  the  CPU  time.  This  problem  does  demonstrate  that  the 
quadratic  method  can  solve  it  in  half  the  number  of  iterations;  therefore,  if  each  analysis 
took  much  longer,  then  a  savings  could  be  achieved.  This  would  the  case  in  practical 
design  problems  when  the  finite  element  model  has  thousands  of  degrees  of  freedom. 
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Figure  5-12:  12  Member  Frame  Reciprocal  Method  Optimization 
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5. 2. 4. 1.2  Reciprocal  Variable  Optimization 

Figure  5-12  shows  the  comparison  of  the  four  reciprocal  methods.  Table  5-18  gives  the 
exact  number  of  iterations  to  converge  to  the  required  tolerance  and  the  CPU  time. 

Table  5-18:  12  Member  Frame  Convergence 


Method 

No.  Of  Iterations  to  Converge 

Time  In  CPU  Sec 

MQA-R 

11 

114.7 

SPA-RS 

14 

77.8 

SPA-RB 

17 

95.4 

SQA-R 

8 

58.0 

The  two  quadratic  methods  solved  this  problem  in  fewer  iterations  than  the  single  point 
first  order  reciprocal  methods;  however,  from  the  graph  one  can  see  that  all  the  methods  are 
very  close.  It  just  takes  the  first  order  methods  a  httle  longer  to  hone  in  on  the  answer 
within  tolerance.  It  is  probably  reasonable  to  conclude  that  the  increased  accuracy  of  the 
quadratic  method  speeds  the  final  convergence. 
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5.2.5  ACOSS  II  (6  Design  Variables) 


5.2.5. 1  Description 

The  Active  Control  of  Space  Structures  (ACOSS)  model  n  was  developed  by  the 
Charles  Stark  Draper  Laboratory  (Canfield,1990:1121).  The  structure  consists  of  two 
subsystems:  the  optical  support  structure  and  the  equipment  section.  The  two  are  connected 
by  springs  at  three  points  to  allow  vibration  isolation.  In  this  problem  the  equipment  section 
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at  the  base  was  disregarded  and  only  the  optical  support  structure,  fixed  at  the  three 
connection  points,  was  considered.  The  finite  element  model  for  this  modified  ACOSS  II 
(with  supporting  cross-members  added)  shown  in  Figure  5-13  has  33  nodes,  18 
concentrated  masses,  and  113  rod  elements  made  of  graphite  epoxy  with  a  Young’s 
Modulus  of  18.5  Mpsi,  weight  density  of  0.055  Ib/in^  and  initial  areas  of  10.0  in^  for  the 
truss  members.  For  the  purposes  of  this  problem  the  structure  was  hnked  to  6  design 
variables  subject  to  a  lower  bound  frequency  of  2.0  Hz,  a  second  frequency  constraint  of 
3.0  Hz,  and  minimum  sizes  of  0.1  in^ 

5. 2. 5. 2  Results 

5. 2. 5. 2.1  Optimization  Problem  Overview 

Table  5-19  summarizes  the  optimization  problem  for  ACOSS  II  hnked  to  six  design 
variables.  Table  5-20  shows  how  the  113  elements  were  linked.  The  approximate  problem 
proved  difficult  to  solve  by  “constr”  in  MATLAB  for  all  the  cases.  Much  trial  and  error 
went  into  finding  the  proper  move  limit  so  that  all  the  cases  could  be  compared  on  an  equal 
basis. 


Table  5-19:  ACOSS  II  (6)  Optimization  Overview 


No.  of  Design  Variables 

6 

No.  of  constraints  (type) 

2(Freq.) 

Convergence  Criteria 

Max.  Change  in  Objective 

5 

Function 

Max.  Constraint  Violation 

.002 

Max.  Change  in  Design  Variable 

.2 

Move  Limits 

.7  with  .9  reduction 
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Min.  Allowable  Design  Variable 

.1 

Max.  Allowable  Design  Variable 

40 

Initial  Objective  Function 

17610 

Initial  Design  Variables 

10,10,10,10,10,10, 

Initial  Max.  Constraint  Violation 

.637 

Optimum  Objective  Function 

13014 

Optimum  Design  Variable  Values 

.22,  13.38,  33.48, 

5.21,  .49,  1.40 

Table  5-20:  ACOSS  II  (6)  Linking 


Design 

Variable 

Elements  linked  to  this  Design  Variable 

1 

1  2  10  11  20  23  33  34  68  69  70  71  94  98  100  104  86  88  106  107  108  109 

2 

4  6  44  46  48  49  52  54  21  22 

3 

3  9  18  19  24  25  28  29  36  37  26  27  30  35 

4 

31  32  12  13  14  15  16  17  45  47  50  51  53  55  56  58  60  62  64  66  57  59  61  63  65 
67  72  73  76  77 

5 

38  39  40  41  42  43  78  80  90  92  79  93  81  82  83  84  85  87  95  97  101  103  99  105 

6 

74  75  89  91  96102  110  111  112  113 

5-31 


NORMALIZED  NORMALIZED 

CONSTRAINT  VIOLATION  OBJECTIVE  FUNCTION 


ACOSS  II  LINKED  TO  6  DV'S 


ACOSS  II  LINKED  TO  6  DV'S 


Figure  5-14:  ACOSS  II  (6)  Direct  Method  Optimization 


5. 2. 5. 2. 2  Direct  Variable  Optimization 

Figure  5-14  shows  the  results  for  the  three  methods  using  direct  design  Table  5-21 
shows  the  exact  number  of  iterations  each  method  took  to  reach  the  convergence  criteria. 
This  chart  also  shows  the  CPU  time  each  method  required  to  converge. 

Table  5-21:  ACOSS  II  (6)  Convergence 


Method 

No.  Of  Iterations  to  Converge 

CPU  Sec 

constr 

42 

1146.4 

MQA-D 

22 

1036.8 

SPA-DS 

25 

819.8 

MQA-D  and  SPA-DS  effectively  solved  the  problem  in  20  iterations,  while  “constr”  took 
considerably  longer  to  get  close.  For  this  problem  the  expense  of  the  quadratic 
approximation  paid  off  with  MQA-D  converging  in  only  22  iterations.  SPA-DB  could  not 
provide  a  solution  for  the  first  iteration. 
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5. 2. 5. 2. 3  Reciprocal  Variable  Optimization 


Figure  5-15  shows  the  comparison  of  the  four  reciprocal  methods.  Table  5-22  gives  the 
exact  number  of  iterations  and  the  CPU  time  to  converge  to  the  required  tolerance. 

Table  5-22:  ACOSS  II  (6)  Convergence 


Method 

No.  Of  Iterations  to  Converge 

CPU  Sec 

MQA-R 

19 

1068.9 

SPA-RS 

19 

864.3 

SPA-RB 

20 

938.3 

SQA-R 

25+ 

908.8 

The  reciprocal  variables  did  just  shghtly  better  than  the  direct  for  this  problem.  It  was 
interesting  that  this  is  first  case  where  both  of  the  quadratics  did  not  show  any 
improvement  over  the  first  order  approximations.  SQA-R  had  some  difficulty  with  this 
problem,  stopping  at  25  iterations,  about  150  lb.  above  the  optimum. 

5.2.6  ACOSS  II  (31  Design  Variables) 

5. 2. 6.1  Description 

The  ACOSS  n  model  here  is  exactly  the  same  as  the  one  above  except  that  the  structure  is 
linked  to  31  design  variables.  Table  5-24  shows  how  the  113  elements  are  linked  to  31. 

5. 2. 6. 2  Results 

5. 2. 6. 2.1  Optimization  Problem  Overview 

Table  5-23  gives  an  overview  of  the  optimization  problem  for  ACOSS  n  linked  to  31 
Design  variables. 
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Table  5-23:  ACOSS  II  (31)  Optimization  Overview 


No.  of  Design  Variables 

31 

No.  of  constraints  (type) 

2(Freq.) 

Convergence  Criteria 

Max.  Change  in  Objective 

5 

Function 

Max.  Constraint  Violation 

.002 

Max.  Change  in  Design  Variable 

.2 

Move  Limits 

.8  with  .98  reduction 
except  for  MQA-R(1,.98) 

Min.  Allowable  Design  Variable 

.1 

Max.  Allowable  Design  Variable 

40 

Initial  Objective  Function 

17610 

Initial  Design  Variables 

10  for  all  31 

Initial  Max.  Constraint  Violation 

.637 

Optimum  Objective  Function 

11705 

Optimum  Design  Variable  Values 

.10, 12.65, 21.08, 34.53, .10, 35. 14, 7.29, 
4.93, 12.37,3.09,.70, 10.98,27.47,32.18, 
.10,4.09,5. 19,.  10,2. 84, 1.78, .25, .1,1.21, 
.28,.5836,.10,.56,.45,.10,.95 
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Table  5-24:  ACOSS  II  (31)  Linking 


Design 

Variable 

Elements  linked  to  this  Design  Variable 

1 

12  10  11 

2 

46 

3 

39 

4 

24  25  18  19 

5 

20  23 

6 

28  29  36  37 

7 

31  32 

8 

12  13  14  15  16  17 

9 

44  46  48  49  52  54 

10 

45  50  47  51  55  53 

11 

38  39  40  41  42  43 

12 

21  22 

13 

26  27 

14 

30  35 

15 

33  34 

16 

56  58  60  62  64  66 

17 

57  59  61  63  65  67 

18 

68  69  70  71 

19 

72  73  76  77 

20 

74  75  89  91 

21 

78  80  90  92 

22 

98  94  100  104 

23 

96  102 

24 

79  93 

25 

81  82  83  84 

26 

86  88 

27 

85  87 

28 

95  97  101  103 

29 

99  105 

30 

106  107  108  109 

31 

no  111  112 
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NORMALIZED  NORMALIZED 

CONSTRAINT  VIOLATION  OBJECTIVE  FUNCTION 


ACOSS  II  LINKED  TO  31  DV'S 


ACOSS  II  LINKED  TO  31  DV'S 


Figure  5-16:  ACOSS  11  (31)  Direct  Method  Optimization 
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5. 2. 6. 2. 2  Direct  Variable  Optimization 

Figure  5-16  shows  the  results  for  the  three  methods  using  direct  design  variables.  Also 
shown  in  Table  5-25  are  the  exact  number  of  iterations  and  the  CPU  time  each  method  took 
to  reach  the  convergence  criteria. 

Table  5-25:  ACOSS  II  (31)  Convergence 


Method 

No.  Of  Iterations  to  Converge 

CPU  Sec 

constr 

200 

5712 

MQA-D 

36 

14896 

SPA-DS 

49+ 

3919.2 

As  with  ACOSS  II  linked  to  six  variables  this  problem  proved  difficult  to  solve  and 
several  trial  and  error  cases  were  needed  to  find  the  proper  move  limit  so  that  the  optimizer 
in  MATLAB  could  solve  the  approximate  problem.  As  with  the  six  variable  case  SPA-DB 
would  get  stuck  on  the  second  iteration.  MQA-D  has  a  clear  advantage  in  solving  this 
problem  when  only  the  number  of  iterations  is  considered.  It  gives  useful  results  in  about 
20  iterations.  However,  to  let  it  proceed  past  this  point  cost  a  substantial  amount  of  CPU 
time.  The  large  amount  of  CPU  time  required  to  solve  the  problem  in  the  quadratic  case  is 
related  to  the  fact  that  the  finite  element  analysis  has  to  solve  a  problem  with  approximately 
90  degrees  of  freedom;  whereas,  the  pseudo-inverse  calculations  involve  finding  496 
unknowns  in  the  Hessian.  The  Hessian  calculation  is  expected  to  be  much  more 
expensive.  However,  if  the  problem  had  31  design  variables  and  many  thousands  of 
degrees  of  freedom  then  estimating  the  Hessian  would  not  have  used  the  majority  of  the 
CPU  time  and  MQA-D  would  have  been  more  economical. 
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Figure  5-17:  ACOSS  II  (31)  Reciprocal  Optimization 
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5. 2. 6. 2. 3  Reciprocal  Variable  Optimization 


Figure  5-17  shows  the  comparison  of  the  four  reciprocal  methods.  Table  5-26  gives  the 
exact  number  of  iterations  and  the  CPU  time  to  converge  to  the  required  tolerance. 

Table  5-26:  ACOSS  II  (31)  Convergence 


Method 

No.  Of  Iterations  to  Converge 

CPU  Sec 

MQA-R 

26* 

8464 

SPA-RS 

28 

2590 

SPA-RB 

86 

N/A 

SQA-R 

52+ 

4635 

*  Used  different  move  limits 


MQA-R  solved  this  problem  in  the  least  iterations  but  it  required  a  different  move  limit  than 
the  rest  for  MATLAB  to  solve  the  approximate  problem.  The  move  limit  began  at  1  and 
was  reduced  by  .98  at  each  iteration.  This  move  limit  would  not  work  for  the  others.  SPA¬ 
RE  came  close  to  the  optimum  in  about  15  iterations  but  it  could  not  converge  until  86 
iterations.  It  just  bounced  around  the  optimum  until  the  move  limit  was  reduced  small 
enough  to  force  it  to  converge.  SQA-R  took  a  very  long  and  conservative  route  to  the 
optimum.  It  was  close  in  about  25  iterations  but  still  did  not  meet  the  convergence  tolerance 
in  52  iterations.  Since  the  move  Umits  are  different  it  is  not  really  fair  to  compare,  but 
MQA-R  does  solve  this  problem  in  fewer  iterations,  even  though  it  used  a  great  deal  of 
CPU  time.  However,  the  two  single  point  reciprocal  methods  are  very  close  for  a  useable 
solution. 
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5.2.7  40  Member  Frame 


Figure  5-18:  40-Member  Frame 


5. 2. 7.1  Description 

The  40-member  plane  frame  shown  in  Figure  5-18  is  made  up  of  I-section  members  with 
32  nodes  with  nodes  3 1  and  32  restrained  from  displacement.  It  was  subjected  to  3  load 
cases  with  288  displacement  constraints  and  120  stress  constraints  shown  in  Table  5-27. 
For  this  stmcture  there  is  a  fixed  relationship  between  the  moment  of  inertia®,  section 
modulus  (S)  and  the  cross  sectional  area  (A).  I=0.20720xA^,  S=0.39300xA^ 
(Kolonay,1987: 109-1 15),  where  the  constant,  0.20720,  has  units  of  in^  and  the  constant, 
0.39300,  has  units  of  in  '.  Material  and  element  properties  are  given  in  Table  5-28  and 
Table  5-29. 


Table  5-27:  40-Member  Frame  Load  Cases 


Direction  of 

oad 

Load  Case 

Node 

x(lb) 

y(ib) 

1 

1 

2 

-24000  1 

3 

4 

-12000  1 

5 

6 

-12000 

7 

-12000 

8 

9 

-12000 

10 

-12000 

11 

12 

-12000 

13 

-12000 

14 

15 

-12000 

16 

-12000 

17 

18 

-12000 

19 

-12000 

20 

21 

-12000 

22 

-12000 

23 

24 

-12000 

25 

-12000 
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Table  5-28:  40-Member  Frame  Material  Properties 

Material _ Steel  _ 

Modulus  of  Elasticity _  29x10°  psi 

Specific  Weight _  .283  Ib/in _ 

Upper  Limit  on  Area _  88.0  in"^ 

Lower  Limit  on  Area _ 2rin _ 

Number  of  Loading  Conditions  3  _ _ 

Displacement  limits(All) _  x(2  in'^),  y(10  in^),z(50.rad) 


Table  5-29:  40-Member  Frame  Element  Properties 


Area 

Section  Modulus 

353.7  in" 

Moment  of  Inertia 

Maximum  Allowable  Stress 

Initial  Weight 

509401b 
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5. 2. 7. 2  Results 


Table  5-30  gives  an  overview  of  the  information  used  in  this  optumzation  problem.  AH 
of  the  methods  converged  to  the  same  optimum. 

Table  5-30:  40-Member  Frame  Optimization  Overview 


No.  of  Design  Variables 

40 

No.  of  constraints  (type) 

288(Disp),  120(Stress) 

Convergence  Criteria 

Max.  Change  in  Objective 

10 

Function 

Max.  Constraint  Violation 

.002 

Max.  Change  in  Design  Variable 

.5 

Move  Limits 

1  with  .9  reduction 

Min.  Allowable  Design  Variable 

.5 

Max.  Allowable  Design  Variable 

88 

Initial  Objective  Function 

50940 

Initial  Design  Variables 

30  for  all 

Initial  Max.  Constraint  Violation 

-.0969 

Optimum  Objective  Function 

34894 

Optimum  Design  Variable  Values 

12.41,12.41,16.97,16.97,19.60,19.66, 

21.72,21.70,23.21,23.24,24.40,24.40, 

25.29,25.38,26.33,26.31,26.30,26.28, 

15.89,15.89,10.60,10.60,12.88,12.92, 

14.75,14.65,16.02,16.00,17.13,17.22, 

18.85,18.83,20.54,20.58,22.32,22.30, 

24.64,24.60,37.86,37.66 
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Figure  5-19:  40-Member  Frame  Direct  Optimization 
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5. 2. 7. 2.1  Direct  Variable  Optimization 

Figure  5-19  shows  the  results  for  the  three  methods  using  direct  design  variables.  Table 
5-31  shows  the  exact  number  of  iterations  each  method  took  to  reach  the  convergence 
criteria.  This  table  also  shows  the  CPU  time  each  method  required  to  converge. 

Table  5-31:  40-Member  Frame  Convergence 


Method 

No.  Of  Iterations  to  Converge 

Time  In  CPU  Sec 

constr 

42 

428 

MQA-D 

9+ 

16378.5 

SPA-DS 

48 

8309.7 

SPA-DB 

43 

1149.9 

MQA-D  and  SPA-DS  were  able  to  reduce  the  objective  function  to  near  optimum  within 
about  five  iterations,  but  they  both  had  trouble  meeting  the  constraint  tolerance.  MQA-D 
was  only  able  to  run  nine  iterations  because  of  the  time  each  iteration  was  taking.  Table  5- 
31  shows  the  amount  of  CPU  time  required  was  much  greater  than  the  others.  Also,  the 
timp,  was  increasing  greatly  for  each  iteration.  This  problem  demonstrates  a  weakness  of 
this  method.  The  amount  of  time  to  calculate  the  approximate  Hessian  for  this  method  has 
become  prohibitively  expensive  as  the  number  of  design  variables  and  previous  points 
increases. 
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Figure  5-20:  40-Member  Frame  Reciprocal  Optimization 
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5. 2. 7. 2. 2  Reciprocal  Variable  Optimization 

Figure  5-20  shows  the  comparison  of  the  four  reciprocal  methods.  Table  5-32  gives  the 
exact  number  of  iterations  to  converge  to  the  required  tolerance  and  the  CPU  time. 

Table  5-32:  40-Member  Frame  Convergence 


Method 

No.  Of  Iterations  to  Converge 

Time  In  CPU  Sec 

MQA-R 

8+ 

11755 

SPA-RS 

10 

742.6 

SPA-RB 

26+ 

1413.4 

SQA-R 

8 

903.3 

As  with  MQA-D  above,  MQA-R  has  become  prohibitively  expensive  for  this  problem. 
It  is  doing  a  very  good  job  for  each  iteration.  In  eight  iterations  it  has  almost  converged  but 
each  iteration  takes  an  unreasonable  amount  of  time  relative  to  the  expense  of  the  exact 
analysis.  This  is  expected,  since  the  finite  element  analysis  is  solving  for  90  degrees  of 
freedom  versus  a  pseudo-inverse  calculation  for  the  820  unknown  terms  in  the 
approximate  Hessian.  SPA-RS  and  SQA-R  both  do  just  as  well  or  better  using  much  less 
CPU  time.  The  reciprocal  methods  for  this  problem  converge  so  quickly  that  very  little 
Hessian  information  is  built  up  for  the  few  previous  points  as  compared  to  the  number  of 
design  variables.  This  problem  confirms  the  presumption  that  MQA  is  probably  only 
useful  for  problems  with  few  design  variables  and  constraints. 
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6.  Conclusion 


6.1  Overview 

Table  6-1  is  a  summary  of  the  number  of  iterations  needed  to  meet  the  convergence 
criteria.  Several  comparisons  can  be  made  from  this  table.  It  should  be  emphasized  that 
the  purpose  of  the  comparison  is  to  look  at  general  trends  in  solving  the  different 
structures.  Other  structures  may  or  may  not  have  similar  results,  but  I  believe  that  these 
results  will  reflect  reasonable  expectations  for  solving  similar  type  problems.  Several 
comparisons  can  be  made,  reciprocal  vs.  direct  variables,  quadratic  vs.  linear  approaches, 
multipoint  vs.  two  point  quadratic  approximations,  and  spherical  vs.  box  move  limits. 


Table  6-1:  Optimization  Summary 


Direct  Variable  Methods 

Reciprocal  Variable  Methods 

constr 

MQA-D 

SPA-DS 

SPA-DB 

MQA-R 

SPA-RS 

SPA-RB 

SQA 

5-Section  Beam 

16 

11 

20 

14 

6 

16 

12 

6 

4-Member  Frame 

47 

38 

116+ 

40 

5 

14 

12 

12 

10  Bar  Truss 

34 

35 

50+ 

73+* 

14 

14 

14 

36+ 

12  Member  Frame 

34 

17 

46 

51 

11 

14 

17 

8 

42 

22 

25 

N/A 

19 

19 

20 

25+ 

ACOSS  11(31) 

200 

36 

49+ 

N/A 

26* 

28 

86 

52+ 

40  Member  Frame 

42 

9+ 

48 

43 

8+ 

10 

26+ 

8 

*  Different  Move  Limits  than  others  for  same  structure 


6.2  Reciprocal  vs.  Direct 

As  expected,  reciprocal  variable  approximations  performed  much  better  in  almost  every 
case.  This  is  typical  with  structural  problems  with  stress  and  displacement  constraints. 
With  the  two  ACOSS  models  the  difference  was  not  as  great,  probably  because  the 
frequency  constraints  tend  to  be  highly  non-linear  and  difficult  to  approximate. 
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6.3  Quadratic  vs.  Linear 


The  focus  of  this  thesis  was  to  test  a  quadratic  approximation,  specifically  MQA.  For 
the  direct  variable  cases  MQA-D  converged  in  the  fewest  number  of  iterations  in  six  out  of 
seven  cases.  For  the  12  member  frame  it  converged  in  less  than  half  the  number  of 
iterations  of  the  other  direct  methods.  MQA-D  does  appear  to  have  a  drawback,  shown 
from  the  cases  where  CPU  times  were  compared.  It  is  computationally  expensive. 
Calculating  the  approximate  Hessian  involves  a  pseudo-inverse  which  solves  for  n(n+l)/2 
elements  in  the  Hessian  matrix.  Where  n  is  the  number  of  design  variables.  Unless  the 
finite  element  analysis  is  of  this  order  or  larger,  then  it  is  doubtful  that  any  benefit  will  be 
obtained  with  regard  to  the  overall  computational  expense  of  solving  the  optimization 
problem.  The  40  Member  Frame  and  ACOSS  hnked  to  31  design  variables  demonstrate 
this  because  for  these  problems  the  Hessian  is  much  more  expensive  to  calculate  than  the 
finite  element  analysis.  However,  if  the  number  of  design  variables  and  constraints  are  not 
too  large  MQA  could  definitely  prove  worthwhile  when  function  evaluations  are  very 
expensive,  since  MQA  for  direct  variables  does  show  promise  in  solving  problems  in  fewer 
iterations.  For  the  reciprocal  variables  MQA  still  did  very  well,  but  as  it  was  already 
known  that  reciprocal  variables  generally  do  a  very  good  job  in  approximating  structural 
problems,  very  little  was  gained  in  most  cases  by  using  the  quadratic  approximation.  It  did 
help  with  accuracy  when  zeroing  in  on  convergence  criteria.  The  same  drawback  of  high 
computational  expense  for  large  numbers  of  variables  found  in  the  MQA  results  in  direct 
variables,  applies  to  MQA  for  reciprocal  variables. 

6.4  Sequential  Two  Point  Quadratic  vs.  Multipoint  Quadratic 

Of  the  cases  tested,  the  multipoint  was  quicker  to  converge  in  three,  slower  in  three,  and 
the  same  in  one  as  compared  to  the  sequential.  From  this  data,  the  only  conclusion  is  that 
using  multiple  points  probably  does  not  increase  the  accuracy  of  the  approximation  greatly 
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over  using  sequential  information  updated  based  on  the  two  most  recent  points.  However, 
overall  the  multipoint  approximation  was  much  more  stable  in  converging  than  the 
sequential.  SQA-R  converged  slowly  for  three  of  the  seven  problems. 

6.5  Spherical  Vs.  Box  Move  Limits 

It  is  not  completely  fair  to  compare  the  two  move  limit  types  without  some  criteria  to 
equate  an  equal  spherical  move  limit  to  some  box  limit.  No  good  way  to  exactly  compare 
these  has  yet  been  developed,  but  from  the  data  obtained  here,  it  appears  that  whether  or 
not  one  type  is  better  than  the  other  is  very  case  dependent.  The  two  ACOSS  problems 
would  not  run  properly  with  box  move  limits  and  the  10  bar  truss  had  to  be  solved  with 
smaller  move  limits  using  box  limits,  because  with  the  same  move  limit  as  the  other  direct 
methods  it  would  diverge.  Generally,  it  appears  that  the  spherical  move  limits  result  in  a 
more  stable  if  not  faster  convergence  for  most  of  the  problems  as  shown  by  the  graphs  in 
chapter  5. 

6.6  Recommendations 

The  next  step  in  testing  MQA  is  to  not  use  every  previous  point  but  to  pick  a  group  of  the 
most  recent  or  closest  points  to  see  if  this  will  increase  the  computationally  efficiency 
without  sacrificing  accuracy.  Also,  if  SQA  strengths  and  weaknesses  can  be  predicted  it 
may  be  very  useful  to  solve  certain  types  of  problems  or  perhaps  it  can  be  modified  to 
perform  well  in  more  cases. 

6.7  The  Bottom  Line 

MQA  does  do  a  good  job  of  approximating  the  functions  to  be  used  in  optimization. 
How  much  better  is  dependent  on  the  particular  problem.  In  aU  cases,  its  convergence  was 
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very  stable  and  it  would  seem  to  be  very  reliable.  However,  whether  or  not  the 
eomputational  expense  of  calculating  the  second  order  approximation  using  MQA  is 
computationally  efficient  will  depend  on  the  relative  size  of  the  analysis  and  the 
optimization  problems.  MQA  is  recommended  only  for  problems  with  few  design 
variables  and  very  expensive  analysis  of  functions  (more  expensive  than  a  singular  value 
decomposition,  which  is  the  majority  of  the  expense  in  calculating  the  pseudo-inverse). 
Also,  it  is  recommended  for  problems  where  appropriate  intermediate  variables,  such  as 
reciprocal  variables  for  structural  problems,  are  not  known. 
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Appendix  A:  Multipoint  Quadratic  Approximation  Matiab  Routine 

function [x, OPTIONS] =mqa (FUN, XO , OPTIONS, VLB, VUB, GRADFUN,APX_TYPE, options) 

%  MQA  Finds  the  constrained  minimiim  of  a  function  of  several  variables 
%  by  Multipoint  Quadratic  Approximations . 

% 

%  X=MQA  (  '  FUN '  ,  XO )  starts  at  XO  and  finds  a  constrained  minimiim  to 

%  the  function  which  is  described  in  FUN  (usually  an  M-file:  FUN.M) . 

%  The  function  'FUN'  should  return  two  arguments:  a  scalar  value  of  the 

%  function  to  be  minimized,  F,  and  a  matrix  of  constraints,  G; 

%  [F,G] =FUN(X) .  F  is  minimized  such  that  G  <  zeros (G) . 

% 

%  X=MQA( ' FUN' ,X, OPTIONS)  allows  a  vector  of  optional  parameters  to 

%  be  defined.  For  more  information  type  HELP  FOPTIONS . 

%  OPTIONS (15)  =  move  limit  reduction  factor  >1  (default=l) 

%  OPTIONS (18)  =  move  limit  ratio  (default=0.5  =  50%) 

% 

%  X=MQA( 'FUN' ,X, OPTIONS, VLB, VUB)  defines  a  set  of  lower  and  upper 

%  bounds  on  the  design  variables,  X,  so  that  the  solution  is  always  in 

%  the  range  VLB  <  X  <  VUB. 

% 

%  X=MQA( 'FUN' ,X, OPTIONS, VLB, VUB, 'GRADFUN' )  allows  a  function 

%  'GRADFUN'  to  be  entered  which  returns  the  partial  derivatives  of  the 

%  function  and  the  constraints  at  X:  [gf,GC]  =  GRADFUN(X). 

%  X=MQA ( ' FUN ' , X , OPTIONS , VLB , VUB , ' GRADFUN ' , APX_TYPE ) 

%  specifies  approximation  type  for  each  function  (obj  &  constraints)  where 

%  1  =  1st  order 

%  2  =  Quadratic 

%  3  =  1st  order  Reciprocal 

%  4  =  Quadratic  Reciprocal 

alpha=l 

%  Written  by  Michael  A. Blaylock  &  Robert  A.  Canfield 

%  Process  input  arguments, 
fcnstr  =  [FUN,  ' (x) ' ] ; 
ndv  =  length (XO); 

defaultOPTIONS  =  [1  0.001  0.001  l.e-4]; 
if  nargin  <  3, 

OPTIONS=def aultOPTIONS ; 

else 

m  =  length (OPTIONS) ; 
for  n=length (defaultOPTIONS) :-l:2, 
if  m<n,  OPTIONS (n)=0;  end; 

if  OPTIONS (n) ==0,  OPTIONS(n)  =  defaultOPTIONS (n) ;  end; 

end 

if  m<l,  OPTIONS (1)  =  1;  end 

end; 

OPTIONS  =  f options (OPTIONS) ; 
if  OPTIONS (14) ==0,  OPTIONS (14) =ndv*10;  end 
if  OPTIONS ( 15 )==0,  OPTIONS ( 15 )=1;  end 
if  OPTIONS ( 18 )==0,  OPTIONS (18) =0.5;  end 
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if  nargin  <  4,  VLB=[];  end 
if  nargin  <  5,  VUB=[];  end 
if  nargin  <  6,  GRADFUN=[];  end 
if  length (GRADFUN) 

grdstr  =  [GRADFUN,  ' (x) ' ] ; 

else 

error (' Finite  Difference  Gradient  not  yet  implemented') 

end 

i f  nargin  <  7 ,  APX_TYPE= [ ] ;  end 

%  Set  options  for  approximate  problem  relative  to  exact,  unless  given, 
if  nargin  <  8, 

options  =  OPTIONS; 

options (1)  =  0; 

options (4)  =  OPTIONS (4) 710, • 

else 

m  =  find(options==0) ; 
options (m)  =  OPTIONS (m) ; 

if  length (options ) <4  |  anY(m==4) ,  options (4)  =  OPTIONS (4) /lO;  end; 

end 

if  OPTIONS (1)==2,  options (1)  =  1;  end; 

%  Error  check  the  lower/upper  bounds  and  initial  guess. 
lenvlb=length(VLB) ;  if  lenvlb,  xlb  =  VLB(:);  end; 
lenvub=length(VUB) ;  if  lenvub,  xub  =  VUB(:);  end; 
xlb (lenvlb+1 :ndv, 1)  =  -Inf * ones (ndv- lenvlb, 1) ; 
xub(lenvub+l:ndv, 1)  =  Inf *ones (ndv-lenvub, 1) ; 
if  lenvlb* lenvub>0 

range  =  1 : min ( lenvlb, lenvub) ; 

if  any (VLB (range) >VUB (range) ) ,  error (' Bounds  Infeasible'),  end 

end 

X  =  XO ( : ) ; 

dx  =  zeros (size (x) ) ; 

if  any(x<xlb)  [  any(x>xub) 

disp (' Initial  X  vector — not  within  bounds — has  been  reset') 

X  =  max(  X,  xlb  ) ; 

X  =  min(  X,  xub  ) ; 

end 

violated  =  abs (OPTIONS (4) ) ; 
active  =  -10*violated; 
movred  =  abs (OPTIONS ( 15 )) ; 
movmax  =  abs (OPTIONS ( 18 )) ; 
movmin  =  0; 
move  =  movmax ; 

if  movred  <  1,  movred  =  1/movred;  end; 

%  Initial  function  evaluations. 

[f,  g(:)]  =  eval ( fcnstr) ; 

[mg,mj]  =  max(g); 
neon  =  length(g); 

if  length(g)  <  1,  error (' Constraints  must  be  defined'),  end 
gapx=NaN*ones (size (g) ) ; 
bound  =  '  ' ; 

HessF  =  zeros (ndv, ndv); 
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apx_order  =  2*ones (l,ncon+l) ; 

apx_order ( 1 : length {APX_TYPE) )  =  APX_TYPE { : ) ' ; 

apx_f  =  apx_order ( 1 ) ; 

apx_g  =  apx_order (2 :ncon+l) ; 

%  First  order  for  first  iter, 
m  =  find(  apx_g==2  |  apx_g==4  ) ; 
apx_g (m)  =  apx_g (m) - 1 ; 

%  Set  up  iteration  count  and  optional  print, 
niter  =  0; 
countfcn  =  1; 
countgrd  =  0 ; 
if  OPTIONS (1)>0 
disp ( ' ' ) 

disp( ' f-COUNT  FUNCTION  MAX{g}  j  gapx_j  Move_liinit 

inax_dx'  )  ; 
end 


% - Main  Loop - 

converged  =  0 ; 
while  converged  ~=  1 

niter  =  niter  +  1; 
fold  =  f; 
xO  =  x; 

feasible_start  =  mg  <  violated; 
if  OPTIONS {1)>0 

[m,n]  =  max(  abs(dx)  ); 

disp ( [sprintf ( ' %5 . Of  %12.6g  %12.6g  %3.0f  %9.3g  %9.3g  %9.3g',.. 
countfcn, f ,mg,mj ,gapx(mj ) ,move,dx(n) )  bound] ) ; 

end 

%  Solve  quadratic  problem. 

[gradf,  gradg]  =  eval (grdstr) ; 
xstore  =  [xstore,x{ : ) ] ; 
gstore  =  [gstore,g( : ) ] ; 
dgstore  =  [ dgs tore, gradg] ; 
if  apx_f==2 

fstore  =  [fstore,f]; 
dfstore  =  [dfstore, gradf ] ; 

HessF  =  hapx{  fstore,  dfstore,  xstore  ) ; 
if  OPTIONS (1)==3,  HessF,  end 
elseif  apx_f==4 

fstore  =  [fstore, f]; 

dfstore  =  [df store, gradf  .*  (-xO. ''2)  ]  ; 

HessF  =  hapx(  fstore,  dfstore,  1. /xstore  ) ; 

end 

for  j=l:ncon 
dg=  [  ] ; 

if  apx_g(j)==2 

Hess  =  hapx(  gstore (j ,  dgstore (:, j :ncon: j + (niter- 
l)*ncon),  xstore, alpha  ); 

HessG(:,j)  =  Hess ( : ) ; 

if  OPTIONS (1)==3,  Hess,  eig(Hess)',  end 
elseif  apx_g(j)==4 
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dgj  =  dgstore( : , j mcon: j+(niter-l) *ncon) ; 
for  k=l:niter 

dgy  =  dgj  (  :  ,  k)  .  *  (-xstore  { : ,  k)  . ''2)  ; 
dg= [dg,dgy] ; 
end 

Hess  =  hapx{  gstore ( j , : ) ' ,  dg,  1 . /xstore, alpha  ); 

HessG(:,j)  =  Hess(:); 

if  OPTIONS (1)==3,  Hess,  eig(Hess)',  end 

end 

end; 

if  OPTIONS (1)==3,  rank (xstore) ,  end 
%  HessF 

%  Hess 

%  disp('eig  Hess') 

%  eig(Hess) 

dx  =  constr ( ' qapx ' ,  zeros (size (xO) ) ,  options,  xlb-x,  xub-x, 

' qapxgrad ' , . . . 

f,  g,  gradf,  gradg,  HessF,  HessG,  xO,  move,  [apx_f  apx_g]  ) ; 

%  dx=sqp( 'qapx' ,  zeros (size (xO )) ,  options,  xlb-x,  xub-x,  'qapxgrad',... 

%  [],[],£/  g,  gradf,  gradg,  HessF,  HessG,  xO,  move,  [apx_f  apx_g]  ) ; 

X  =  xO  +  dx; 

[fapx,gapx]  =  qapx(  dx,  f,  g,  gradf,  gradg,  HessF,  HessG,  xO,  move,... 

[apx_f  apx_g]  ) ; 
countgrd  =  countgrd  +  1; 
feaslp  =  max(gapx)  <=  options(4); 
if  gapx(ncon+l)  >  active 
bound  =  '  bound ' ; 

else 

bound  =  '  ' ; 

end 

%  Evaluate  functions  at  new  point. 

[f,  g(:)]  =  eval (fcnstr) ; 

[mg,mj]  =  max(g); 
countfcn  =  countfcn  +  1; 

%  Increase  move  limit  if  no  constraints  encountered  and  design  improved, 
if  feasible_start  &  f<fold  &  mg  <  5*active 

if  niter>l,  move  =  move  +  movmax/movred;  end; 

%  Reduce  move  limit  if  linear  approximation  was  feasible, 
else 

move  =  movmin  +  (move-movmin) /movred; 

end 

%  Check  convergence . 

if  max(abs (dx) )<OPTIONS(2)  &  abs (f-fold) <OPTIONS (3 )  &  mg<OPTIONS ( 4 ) 
if  OPTIONS (1)>0 

disp ( [sprintf ( ' %5 . Of  %12 . 6g  %12.6g 
%3 . Of ' , countfcn, f ,mg,mj ) ] ) ; 

disp (' Optimization  Terminated  Successfully') 
Active_Constraints  =  f ind (g>active) 

end 
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OPTIONS (14)  =  niter; 
converged  =  1 ; 
elseif  niter  >=  OPTIONS (14) 

disp (' Maximum  niimber  of  iterations  exceeded  in  MQA' ) 
disp (' increase  OPTIONS (14) ' ) 
converged  =  1 ; 

else 

apx_g  =  apx_order (2 :ncon+l) ;  %  Restore  approximation  type  after 

first  iter, 
end 

end 

%  Return  iteration  history  info  in  OPTIONS. 

OPTIONS (10)  =  countfcn; 

OPTIONS (11)  =  countgrd; 
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Appendix  B:  Routine  For  Calculating  the  Approximate  Hessian 

function  Hess  =  hapx(  f,  g,  x,  alpha,  print  ) 

% 

%  Quadratic  Approximation:  find  approximation  to  Hessian 
%  based  on  current  &  previous  information 

% 

%  INPUTS 
% 

%  f  =  function  values  col\imn  vector 
%  g  =  gradients  (each  column  is  a  gradient  of  f) 

%  X  =  design  points  (each  column  is  a  point  where  f  &  g  evaluated) 

%  alpha  =  restriction  on  negative  terms  0  (no  neg)  to  1  (allows  full  neg) 

%  print  =  optional  print  control 
% 

%  OUTPUT 

%  Written  by  Robert  A.  Canfield 

%  Hess .  Hessian  matrix  approximation  (n  by  n) 

% 

[ndv,npt]  =  size(  x  ) ; 

if  npt  <  2  1  ndv  <  1,  Hess=zeros (length (x) , length (x) ) ;  return;  end; 

% 

%  Local  variables 


% 

%  ndv=n. . . .  Number  of  Design  Variables 

%  npt=k. . . .  Number  of  Points  where  f  &  g  evaluated 

%  neg .  Number  of  Equality  constraints 

%  N .  Number  of  independent  terms  in  Hessian  matrix 

%  p .  Set  of  design  points  where  equality  constraints  met 

%  A .  LHS  coefficient  matrix  of  Hessian  vector 

%  b .  RHS  of  equality  constraints  for  function  values 

%  xO .  Current  (k_th)  design  point 

%  dx .  Change  between  previous  and  current  design  vectors 

%  diagl,n. .  term  in  h  where  each  upper  diagonal  of  H  begins /ends 

%  beta . weighting  factor  for  each  previous  design  point 

% 


%  Initialize  local  variables 
% 

if  nargin<4 
alpha  =1; 
end 
alpha; 

nprev  =  npt-1; 
neq  =  nprev; 

%neq  =  min(  nprev,  ndv  ) ; 

N  =  ndv* (ndv+1) /2 ; 

P  =  1 : nprev ; 

k  =  npt; 

xO  =  x( : ,k) ; 

dx  =  x(:,P)  -  diag(xO) *ones (ndv,neq) ; 
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b  =  f(P)  -  f (k) *ones (neq, 1)  -  dx'*g{:,k); 
b  =  max (b, alpha *b) ; 

% 

%  Create  index  vectors  to  convert  H  from  matrix  to  vector  and  back. 
%  (Upper  diagonals  of  H  are  stored  in  a  vector  h) . 

% 

ndiag  =  ndv; 

diagl  =  ones (1, ndiag) ; 

diagn  =  zeros (1, ndiag) ; 

diagn(l)  =  ndv; 

for  n=2: ndiag, 

diagl (n)  =  diagn (n-1)  +  1; 
diagn(n)  =  diagl (n)  +  ndv-n; 

end; 

% 

%  Transformation  matrix  T  converts  n(n+l)/2  symmetric  h  vector 
%  to  vector  of  rows  of  full  H. 

% 

Trow  =  l:ndv''2; 

Tool  =  zeros  (l,ndv''2)  ; 
t=0; 

for  row=l:ndv, 
for  col=l:ndv, 
t  =  t+1; 

u  =  abs (col-row) ; 

Tcol{t)  =  diagl (u+l)-l  +  max (row, col) -u; 

end; 

end; 

T  =  sparse (Trow, Tool , 1) ; 

% 

%  Weighting  factors  for  each  design  point. 

% 

beta  =  zeros (l,nprev) ; 

for  p=l:nprev,  beta(p)  =  norm(  x(:,p)  -  xO  )''2;  end; 
beta  =  max(  beta  )  ./  beta; 

% 

%  Form  A:  linear  coefficients  of  H  in  constraints,  A  h  =  b. 

%  Each  row  of  A  corresponds  to  one  of  the  P  constraints. 

%  Columns  of  A  are  coefficients  of  H  in  1/2  dx'*H*dx 

%  arranged  by  diagonal  and  then  codiagonal  terms . 

% 

A  =  zeros (neq, N) ; 
for  row=l:neq; 

dx  =  x(:,P(row))  -  xO; 

A  (row,  1:  ndv)  =  (dx.''2  /  2)'; 
for  n  =  2: ndiag, 

col  =  diagl (n) : diagn (n) ; 

I  =  l;l+ndv-n; 

J  =  n:ndv; 

A (row, col)  =  dx(I) '  .*  dx(J)'; 

end; 

end; 

% 

%  Error  check  against  quadratic  equations 
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for  m=l:neg, 
p  =  P  (m)  ; 
dx  =  x( : ,p)  -xO; 
dg  =  g( : ,k)  -  g( : ,p) ; 

% 

%  If  slopes  are  same,  function  values  must  be  consistent. 

% 

%  slopeO  =  g(:,k)'*dx; 

%  slopep  =  g(:,p)'*dx; 

%  if  slopeO*slopep  >=  0  &  abs (slopeO* (f (p)  -  f{k)))  <=  eps, 

%  slopeO,  slopep,  f(k),  f (p) 

%  error (' Bad  slopes') 

%  end 

end;  %  end  for  each  previous  point 

% 

%  Weighted  equality  constraints  for  function  values  and  gradients. 
% 

dx  =  (  x{:,P)  -  xO*ones (l,nprev)  )  *  diag(beta) ; 
dg  =  diag(beta) * (g{ : , P)  -  g ( : , k) *ones (l,nprev) ) ' ; 

%  Magnify  equality  constraints  for  function  values  vs.  gradients, 
magnify  =  max (max (abs (dg) ) )  /  max ( [min (abs (b) ) ,  10^ (-8)]); 

magnify  =  max(  magnify,  1  ) ; 

AA  =  magnify  *  diag (beta) *A; 
bb  =  magnify  *  diag (beta) *b; 

%bb  =  [bb;  dg ( : ) ] ; 

for  p=l:max(P) ; 

if  dg(p, ; ) *dx( : ,p)<0; 

dg(p, : ) =alpha*dg(p, : ) ; 

end 

end 

bb  =  [bb;  dg( : ) ] ; 
for  n=l:ndv, 

DelT  =  zeros(nprev,ndv^2); 

DelT( ; , (n-1) *ndv+l :n*ndv)  =  dx'; 

AA  =  [AA;  DelT*T] ; 

end; 

h  =  pinv(  AA  )  *  bb; 

% 

%  Convert  h  to  H. 

% 

Hess  =  diag(  h(l:ndv)  ); 
for  n=l:ndv-l, 

codiag  =  h(  diagl (n+1) ;diagn(n+l)  ); 

Hess  =  Hess  +  diag (codiag, n)  +  diag (codiag, -n) ; 
end; 

% 

%  Optionally  check  final  residual  and  constraints. 

% 

if  nargin<5,  print=0;  end; 
if  print>0, 
beta 
magnify 
wtdres  =  0 ; 
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constant  =  0; 
for  in=l:neq, 
p  =  P  (m)  ; 
dx  =  X  ( :  ,  p )  -xO; 
dg  =  g( : ,k)  -  g( : ,p) ; 

R  =  Hess*dx  +  dg 

wtdres  =  wtdres  +  beta (p) * (R' *R) /2 ; 
constant  =  constant  +  beta (p) * (dg' *dg) 

end; 

A*h  -  b 
wtdres 

%residual  =  c ' *h  +  h'*Q*h/2  +  constant/2 
end; 
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Appendix  C:  Routine  for  Calculating  the  approximate  Functions 

function  [f,  g]  =  qapx(  dx,  fO,  gO,  gradf,  gradg,  Hf,  Hg,  xO,  move,  apxtype  ) 
% 

%  qapx  Quadratic  Approximations  to  f  and  g 

% 

%- Inputs 
% 

%  dx .  Change  in  independent  variable  vector 

%  fO .  Objective  function  value  at  base  point  dx=0 

%  gO .  Constraint  vector  at  base  point  dx=0 

%  gradf....  Gradient  of  objective  w.r.t.  X 

%  gradg....  Gradients  of  constraints  w.r.t.  X 

%  Hf .  Hessian  of  objective  w.r.t.  X  or  Y  (apxtype=2  or  4,  respectively) 

%  Hg .  Hessian  of  constraints  (each  coliomn  is  vectorized  matrix) 

%  xO .  Original  X  vector 

%  move .  Spherical  Move  limit  (not  applied  if  move=0) 

%  apxtype..  Approximation  Type  (l=linear,  2=quadratic,  3=reciprocal , 

4=rec . quad. ) 

% 

%-Ouputs 

% 

%  f .  Approximate  function  value 

%  g .  Approximate  constraints 

%  (and  spherical  move  limit  constraint  in  last  entry  if  move>0) 

%-Local  Variables 
% 

%  ndv .  Number  of  Design  Variables 

%  neon .  Number  of  Constraints 

%  dr .  Equivalent  perturbation  for  1st  order  Reciprocal  approximation 

%  dy .  Change  in  intermediate  (reciprocal)  variable 

%  Hgj .  Hessian  approximation  for  j_th  constraint 

%  scale....  X  relative  scale  factors  for  norm  constraint 
% 

%  Written  by:  Michael  A.  Blaylock  and  Robert  A.  Canfield 
%  Created:  8/12/94 

%  Modified:  2/  4/95 

% 

%  Modifications 

%  2/4/95  -  Handle  nonpositive  move  limit, 
dx  =  dx  (  :  )  ; 

[ndv, neon]  =  size(  gradg  ) ; 
if  nargin<9,  move=l;  end; 
if  nargin<10 

apxtype=2 *ones ( 1 , ncon+1 )  ; 
elseif  any (apxtype>2) 

dr  =  dx. *x0 . / (xO+dx) ; 
dy  =  l./(xO+dx)  -  l./xO; 

end 
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%  Objective  Approximation, 
if  apxtype (1) ==1 

f  =  fO  +  gradf { : ) ' *dx; 
elseif  apxtype (1) ==2 

f  =  fO  +  gradf (:)'*dx  +  dx ' *Hf *dx/2 ; 
elseif  apxtype (1) ==3 

f  =  fO  +  gradf ( : ) ' *dr; 
elseif  apxtype (1) ==4 

f  =  fO  +  gradf {;)' *dr  +  dy'*Hf*dy/2; 

end 

%  Loop  for  each  constraint  to  unpack  Hessian, 
for  j=l:ncon 

if  apxtype (j +1) ==1 

g(j)  =  gO(j)  +  gradg( : , j ) ■ *dx; 
elseif  apxtype (j+1) ==2 

Hgj  =  reshape {  Hg{:,j),  ndv,  ndv  ); 
g(j)  =  gO{j)  +  gradg { : , j ) ' *dx  +  dx' *Hgj *dx/2 ; 
elseif  apxtype (j +1) ==3 

g(j)  =  gO(j)  +  gradg(:, j) '*dr; 
elseif  apxtype (j+1) ==4 

Hgj  =  reshape (  Hg ( ; , j ) ,  ndv,  ndv  ) ; 

g(j)  =  gO(j)  +  gradg {:, j) '*dr  +  dy ' *Hgj *dy/2 ; 

end 

end 

%  Add  Euclidean  norm  move  limit 
if  move>0 

scale  =  xO; 
small  =  sqrt(eps); 
for  i=l:ndv 

if  scale(i)<small,scale(i)=l;  end; 

end; 

g(ncon+l)  =  (  norm  (dx. /scale)  /  move  ) ''2  -  1; 

end 
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Appendix  D:  Routine  for  Calculating  the  Approximate  Gradients 


function  [df,  dg]  =  qapxgrad(  dx,  fO,  gO,  gradf,  gradg,  Hf,  Hg,  xO,  move. 


apxtype  ) 

%qapxgrad  Quadratic  Approximations  to  f  and  g 
% 

%  dx .  Change  in  independent  variable  vector 

%  fO .  not  used 

%  gO .  not  used  (required  for  consistency  with  calling  sequence  of  qapx) 

%  gradf....  Gradient  of  objective  w.r.t.  X 
%  gradg....  Gradients  of  constraints  w.r.t.  X 

%  Hf .  Hessian  of  objective  w.r.t.  X  or  Y  (apxtype=2  or  4,  respectively) 

%  Hg .  Hessian  of  constraints  (each  column  is  vectorized  matrix) 

%  xO .  Original  X  vector 

%  move .  Spherical  Move  limit 

%  apxtype..  Approximation  Type  (l=linear,  2=quadratic,  3=reciprocal , 


4=rec . quad. ) 

%  Local  Variables 


% 

%  ndv .  Number  of  Design  Variables 

%  neon . Number  of  Constraints 

%  mag .  Magnification  for  1st  order  Reciprocal  approximation  gradients 

%  dy .  Change  in  intermediate  (reciprocal)  variable 

%  dydx .  Chain  rule  for  reciprocal  to  direct  variables 

%  Hgj .  Hessian  approximation  for  j_th  constraint 

%  scale....  X  relative  scale  factors  for  norm  constraint 
% 

%  Written  by:  Michael  A.  Blaylock  and  Robert  A.  Canfield 
%  Created:  8/12/94 

%  Modified:  2/  6/95 


%  Modifications 

%  2/6/95  -  Handle  nonpositive  move  limit. 


% — BEGIN 
% 

%  Initialize  variables. 

[ndv, neon]  =  size(  gradg  )  ; 
dx  =  dx  ( : )  ,- 

if  nargin<9,  move=l;  end; 
if  nargin<10 

apxtype=2*ones (l,ncon+l)  ; 
elseif  any (apxtype>2 ) 

mag  =  (xO  .  /  (xO+dx)  )  .  ^^2  ; 
dy  =  l./(xO+dx)  -  l./xO; 
dydx  =  -  (xO+dx)  . (-2)  ; 

end 

if  move>0 

dg  =  zeros (ndv, ncon+1) ; 

else 

dg  =  zeros (ndv, neon); 
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end 


%  Objective  Approximation  Gradient, 
if  apxtype (1) ==1 

df  =  gradf ( : ) ; 
elseif  apxtype (1) ==2 

df  =  gradf {: )  +  Hf*dx; 
elseif  apxtype (1) ==3 

df  =  gradf (:). *mag; 
elseif  apxtype (1) ==4 

df  =  gradf (:). *mag  +  (Hf *dy) . *dydx; 

end 

%  Loop  for  each  constraint  to  unpack  Hessian, 
for  j=l:ncon 

if  apxtype (j+1) ==1 

dg ( : , j )  =  gradg ( : , j ) ; 
elseif  apxtype ( j+1) ==2 

Hgj  =  reshape (  Hg(:,j),  ndv,  ndv  ); 
dg(:,j)  =  gradg (:,j)  +  Hgj*dx; 
elseif  apxtype { j+1) ==3 

dg { : , j )  =  gradg ( : , j ) . *mag; 
elseif  apxtype (j+1) ==4 

Hgj  =  reshape (  Hg(:,j),  ndv,  ndv  ); 
dg(:,j)  =  gradg (:, j ). *mag  +  (Hgj *dy) . *dydx; 

end 

end 

%  Euclidean  norm  move  limit  gradient, 
if  move>0 

scale  =  xO; 
small  =  sgrt(eps); 
for  i=l:ndv 

if  scale (i)<small, scale (i)=l;  end; 

end; 

dg(:,ncon+l)  =  2  * (dx . /scale . ^2 ) / (move^2 ) ; 

end 
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